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Abstract 


This  report  discusses  finite  element  modeling  of  progressive  damage  in  the  adhesive 
bond  layers  of  actuated  plates  and  investigates  the  reduction  in  actuation  capacity  caused 
by  the  damaged  bond  layers.  The  primary  challenge  posed  by  this  class  of  problems 
stems  from  the  vast  range  of  geometric  scales  that  are  represented,  with  the  thickness  of 
the  adhesive  layer  representing  the  smallest  scale,  the  overall  thickness  of  the  actuated 
plate  representing  the  intermediate  scale  and  the  in-plane  dimensions  of  the  plate 
representing  the  largest  scale.  In  multiscale  problems,  the  overall  efficiency  of  the 
numerical  methodology  is  of  paramount  importance,  thus  model  development  is  guided 
by  the  need  to  obtain  a  sufficiently  accurate  solution  at  an  acceptable  computational 
expense.  In  this  study,  this  goal  is  achieved  by  through  the  use  of  a  hierarchical, 
displacement-based,  2-D  finite  element  model  that  includes  the  first-order  shear 
deformation  model  (FSD),  type-I  layerwise  models  (LW1)  and  type-II  layerwise  models 
(LW2)  as  special  cases.  Both  the  LW1  layerwise  model  and  the  more  familiar  FSD 
model  use  a  reduced  constitutive  matrix  that  is  based  on  the  assumption  of  zero 
transverse  normal  stress;  however,  the  LW1  model  includes  discrete  layer  transverse 
shear  effects  via  in-plane  displacement  components  that  are  C°  continuous  with  respect  to 
the  thickness  coordinate.  The  LW2  layerwise  model  utilizes  a  full  3-D  constitutive  matrix 
and  includes  both  discrete  layer  transverse  shear  effects  and  discrete  layer  transverse 
normal  effects  by  expanding  all  three  displacement  components  as  C°  continuous 
functions  of  the  thickness  coordinate.  The  hierarchical  finite  element  model  incorporates 
a  3-D  continuum  damage  mechanics  model  that  predicts  local  orthotropic  damage 
evolution  and  local  stiffness  reduction  at  the  geometric  scale  represented  by  the 
individual  material  ply  or,  in  the  case  of  layerwise  models,  by  the  individual  numerical 
layer.  The  results  clearly  demonstrate  that  the  resulting  model  can  efficiently  simulate 
progressive  damage  in  the  adhesive  layers.  For  rectangular  actuator  patches,  the  adhesive 
damage  is  highest  near  the  comers  of  the  actuator  and  is  driven  primarily  by  local 
concentrations  in  the  transverse  normal  and  transverse  shear  stresses.  In  contrast  to 
previous  studies  that  have  shown  that  the  inclusion  of  discrete  layer  transverse  normal 
stress  does  not  significantly  influence  the  predicted  global  deformations,  the  present 
study  shows  that  the  transverse  normal  stress  has  a  very  significant  effect  in  the  initiation 
and  progression  of  localized  damage  in  the  adhesive  layers. 

Keywords:  progressive  damage,  adhesive  bonding,  debonding,  actuated  plates,  finite 
elements 
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1.  Introduction  and  Background 

As  evidenced  by  the  large  number  of  works  cited  in  the  recent  reviews  by  Benjeddou 
(2000)  and  Chopra  (2002),  there  have  been  extensive  efforts  to  develop  analytical  models 
and  finite  element  models  of  laminated  composite  structural  components  (i.e.  beams,  plates, 
and  shells)  that  contain  surface- mounted  or  embedded  piezoelectric  actuators  and/or  sensors. 
On  a  purely  mechanical  basis,  these  modeling  efforts  can  be  broadly  classified  into  two 
basic  categories  according  to  the  kinematics  assumed  by  each  model. 

The  first  category  of  models,  known  as  ‘equivalent  single  layer’  models  (or  ESL 
models),  are  distinguished  by  the  use  of  a  displacement  field  that  exhibits  C1  continuity  with 
respect  to  the  laminate  thickness  coordinate,  i.e.,  the  displacement  components  and  their 
thickness  derivatives  are  continuous  through  the  entire  laminate  thickness  dimension.  The 
primary  advantage  of  the  ESL  models  is  their  high  level  of  computational  efficiency  which 
stems  from  the  fact  that  very  few  functions  need  to  be  determined  in  order  to  completely 
define  the  displacement  field.  However,  despite  the  fact  that  ESL  model  solutions  are 
computationally  expedient,  the  basic  assumption  of  a  C1  thickness  continuity  of  the 
displacement  field  is  overly  restrictive  and,  as  discussed  by  Robbins  and  Chopra  (2006), 
poses  two  major  difficulties  in  the  modeling  of  actuated  plates.  First,  the  continuous 
thickness  derivatives  of  the  displacement  components  prevent  the  ESL  models  from 
delivering  accurate  ply- level  stress  and  strain  fields  for  laminates  where  adjacent  material 
layers  differ  significantly  in  their  mechanical  properties  and  thicknesses.  Second,  the  ESL 
models  cannot  accurately  represent  the  diversion  of  actuation  energy  into  localized 
transverse  shear  deformation  near  the  edges  of  the  actuators;  consequently,  the  ESL  models 
consistently  over-predict  the  global  deformation  of  the  actuated  plate  caused  by  the 
actuators. 

Early  efforts  at  developing  ESL  models  of  actuated  plates  were  essentially  adaptations 
of  the  classical  laminate  theory  (or  CLT)  that  is  based  on  the  Kirchhoff  assumption  that 
transverse  normal  material  fibers  remain  straight  and  normal  to  the  curved  midplane  of  the 
structural  component.  Within  the  CLT  framework,  piezoelectrically  actuated  beam  models 
were  developed  by  Crawley  and  deLuis  (1987),  while  early  actuated  plate  and  actuated  shell 
models  were  developed  by  Lee  (1990)  and  Tzou  and  Gadre  (1989)  respectively.  For 
slightly  thicker  laminates,  the  first  order  shear  deformation  theory  (FSD)  provides 
macroscopic  ESL  kinematics  that  are  more  appropriate  since  it  permits  a  rudimentary  gross 
shear  deformation  that  is  assumed  constant  with  respect  to  the  thickness  coordinate.  FSD 
finite  element  models  of  piezoelectrically  actuated  beams,  plates,  and  shells  were  developed 
by  Robbins  and  Reddy  (1991),  Lin  et  al.  (1996),  and  Miller  and  Abramovich  (1995). 

Higher  order  ESL  models  with  full  thermo -electro- mechanical  coupling  were  developed  for 
quasi- static  analysis  by  Chattopadhyay  et  al.  (1999)  and  for  dynamic  analysis  by  Zhou  et  al. 
(2000).  In  both  of  these  studies,  the  in- plane  displacement  components  were  assumed  to  be 
cubic  functions  of  the  thickness  coordinate,  while  transverse  normal  effects  were  neglected 
via  the  assumption  of  zero  transverse  normal  stress. 

The  second  category  of  models,  known  as  ‘discrete  layer’  models  (or  layerwise  models), 
are  distinguished  by  the  use  of  a  displacement  field  that  exhibits  only  C°  continuity  with 
respect  to  the  laminate  thickness  coordinate,  i.e.,  the  displacement  components  themselves 
are  continuous  through  the  entire  laminate  thickness  dimension,  but  their  thickness 
derivatives  are  permitted  to  exhibit  one  or  more  discontinuities  commensurate  with  the  level 
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of  transverse  discretization  employed.  In  practical  terms,  layerwise  models  discretize  the 
laminate  thickness  dimension  into  a  contiguous  set  of  numerical  layers.  The  displacement 
field  is  then  separately  defined  within  each  numerical  layer  in  such  a  way  that  the 
displacement  components  maintain  continuity  across  interlayer  boundaries;  however  their 
thickness  derivatives  are  not  required  to  be  continuous  across  the  interlayer  boundaries. 
Provided  that  the  number  of  numerical  layers  is  greater  than  or  equal  to  the  number  of 
distinct  material  layers,  the  layerwise  displacement  field  is  capable  of  correctly  representing 
the  kinking  and  warping  of  transverse  normal  fibers  that  is  commonly  observed  in 
multilayer  laminates.  This  particular  modeling  capability  becomes  especially  important  in 
actuated  plates  due  to  the  presence  of  adhesive  bond  layers  that  are  typically  much  more 
compliant  than  the  structural  substrate  material  or  the  actuator  material. 

It  should  be  emphasized  that  full  3-D  finite  element  models  are  classified  as  discrete 
layer  models  provided  that  more  than  one  3-D  element  is  used  to  discretize  the  laminate 
thickness  dimension.  If,  on  the  other  hand,  the  entire  thickness  dimension  of  the  laminate  is 
encompassed  in  a  single  layer  of  3-D  finite  elements,  then  the  resulting  model  is  effectively 
an  ESL  model. 

Reddy  (1987,  1989)  generalized  the  layerwise  concept  by  expanding  the  laminate 
displacement  field  through  the  thickness  in  terms  of  a  1-D  Lagrangian  finite  element 
interpolation  where  both  the  number  of  1-D  elements  and  the  polynomial  order  of  the  1-D 
interpolation  functions  are  arbitrary.  This  data  structure  provides  a  very  convenient 
framework  from  which  to  develop  layerwise  finite  element  models  of  beams,  plates  and 
shells.  A  layerwise  model  for  piezoelectrically  actuated  beams  was  first  reported  by 
Robbins  and  Reddy  (1991)  who  used  an  induced  strain  approach  to  approximate  the 
piezoelectric  effect  (i.e.  similar  to  imposing  a  temperature  change  in  the  piezoelectric 
material  only).  Saravanos  and  Heyliger  (1995)  developed  a  layerwise  beam  model  with 
fully  coupled  electrical  and  displacement  fields.  Both  these  studies  clearly  demonstrated 
that  fundamental  differences  exist  between  the  predicted  solutions  from  layerwise  models 
and  ESL  models.  A  layerwise  composite  plate  model  was  reported  by  Robbins  and  Reddy 
(1993),  once  again  using  an  induced  strain  approach  to  approximate  the  piezoelectric  effect. 
They  demonstrated  that  the  resulting  layerwise  plate  model  produced  laminate  solutions  that 
were  equivalent  to  3-D  finite  element  solutions  provided  that  comparable  levels  of 
discretization  were  used.  Heyliger,  et  al.  (1994)  and  Saravanos,  et  al.  (1997)  developed 
layerwise  plate  models  with  full  electro- mechanical  coupling.  Lee  and  Saravanos  (1997) 
extended  the  layerwise  plate  model  to  include  thermo -electro- mechanical  coupling,  while 
Saravanos  (1997)  developed  a  layerwise  shell  model  with  electro- mechanical  coupling.  In 
addition  to  layerwise  finite  element  solutions,  exact  3-D  elasticity  solutions  have  been 
obtained  for  relatively  simple  problems  by  Heyliger  (1997)  and  Gopinathan,  et  al.  (2000). 

The  material  heterogeneity  of  multilayer  actuated  structural  components  inevitably 
gives  rise  to  stress  concentrations  which  typically  occur  at  the  intersection  of  material 
interfaces  and  free  edges.  These  stress  concentrations  primarily  involve  the  transverse 
shear  stress  and  transverse  normal  stress  and  are  most  pronounced  near  the  edges  of  the 
thin  adhesive  layer  where  they  serve  to  transfer  the  load  from  the  actuator  to  the 
underlying  structural  substrate,  or  vice-versa  in  the  case  of  a  bonded  sensor  (Robbins  and 
Reddy,  1996;  Vel  and  Batra,  2000;  Zhang  et  al.,  2003).  These  locally  elevated 
transverse  stresses  are  particularly  important  since  they  are  likely  to  cause  damage  to  the 
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adhesive  material  which  could  become  severe  enough  to  compromise  the  overall 
functionality  of  the  actuator  or  sensor. 

Despite  the  obvious  importance  of  this  particular  issue,  there  have  been  no  reported 
efforts  to  predict  the  evolution  of  damage  within  the  adhesive  bond  layers  of  actuated 
structural  components.  However,  there  has  been  considerable  research  effort  in 
developing  numerical  methods  to  predict  damage  accumulation  in  fiber-reinforced 
composite  laminates,  and  in  the  present  study,  this  body  of  knowledge  is  leveraged  in 
developing  a  numerical  methodology  to  predict  damage  evolution  in  actuated  structural 
components. 

Most  of  the  past  efforts  to  model  the  evolution  of  distributed  microscopic  damage  in 
fiber- reinforced  composite  laminates  have  relied  on  the  use  of  continuum  damage 
mechanics  (or  CDM).  CDM  is  distinguished  by  the  introduction  of  an  internal  tensor 
field  (i.e.  the  damage  tensor)  that  describes  the  distribution,  density  and  orientation  of 
microcracks  (Krajcinovic,  1996;  Lemaitre,  1996;  Skrzypek  and  Ganczarski,  1999).  The 
evolution  of  the  damage  tensor  and  the  resulting  reduction  of  material  stiffness  are 
described  within  the  framework  of  irreversible  thermodynamics. 

Early  efforts  to  apply  CDM  to  composite  laminates  (e.g.,  Talreja,  1985;  Lene,  1986; 
Allen  et  al.,  1987a,b;  and  Lee  et  al.,  1989)  focused  mainly  on  the  development  of  a 
suitable  form  for  the  damage  variable  and  expressing  the  stiffness  degradation  in  terms  of 
the  chosen  damage  variable.  Talreja  (1985)  chose  a  set  of  vectors  to  represent 
microcrack  densities  on  various  directed  planes  in  the  pre- homogenized  composite 
material,  while  Lene  (1986)  used  a  scalar  quantity  to  describe  the  degree  of  fiber/matrix 
debonding  in  the  heterogeneous  microstructure  prior  to  homogenization.  In  both  of  these 
works,  the  damage  evolution  equations  received  minimal  or  no  treatment.  Allen  et  al. 
(1987a,b)  and  Lee  et  al.  (1989)  developed  a  progressive  damage  model  of  laminated 
composites  using  damage  dependent  constitutive  relations  for  the  homogenized 
composite  material  at  the  ply  level.  Distributed  microscopic  damage  was  quantified 
using  an  internal  variable  (2nd  order  tensor)  that  represented  additional  strain  that  results 
from  the  presence  of  damage;  however,  the  damage  evolution  law  was  empirically  based, 
instead  of  being  derived  from  a  free  energy  function.  The  model  accounted  for  matrix 
cracking  and  fiber  fracture  for  laminates  subjected  to  monotonic  and  cyclical  uniaxial 
loading. 

Ladeveze  and  co-workers  (e.g.,  Ladeveze  and  Dantec,  1992;  Allix  and  Ladeveze, 
1989,1992;  Ladeveze  et  al.,  2000)  developed  a  mesomechanical  model  where  damage  is 
predicted  independently  for  each  homogenized  composite  ply  and  each  interface  that 
separates  adjacent  plies.  The  in-plane  damage  modes  are  assumed  to  be  uniformly 
distributed  through  the  thickness  of  each  homogenized  composite  ply  and  are  described 
by  three  internal  variables  that  define  the  reduction  in  the  in-plane  stiffnesses  En,  E22  and 
G12.  The  transverse  damage  mode  is  restricted  to  the  2-D  interfaces  separating  adjacent 
plies  and  is  described  by  an  additional  internal  variable  that  effectively  defines  the 
reduction  in  transverse  normal  and  transverse  shear  stiffnesses  E33,  G13,  G23.  The  damage 
variables  are  assumed  to  have  simple  evolution  laws  that  are  derived  from  a  free  energy 
function. 

Voyiadjis  and  coworkers  (e.g.,  Voyiadjis  and  Kattan,  1993,  1999;  Voyiadjis  and  Park, 
1999;  Voyiadjis  and  Delikas,  2000)  developed  a  3-D  model  for  coupled  progressive 
damage  and  plasticity  based  on  the  use  of  a  symmetric  second  order  damage  tensor 
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whose  the  eigenvectors  represent  the  principal  directions  of  damage  and  whose 
eigenvalues  represent  the  density  of  distributed  microcracks  that  are  normal  to  the 
respective  eigenvectors.  Damage  is  predicted  separately  for  the  fiber  and  matrix 
constituents  followed  by  homogenization  of  the  damaged  microstructure.  Although 
rigorously  formulated  using  a  three- step  split  operator  algorithm  to  separate  the 
deformation  into  elastic,  plastic  and  damage  components,  the  resulting  model  appears  to 
be  quite  computationally  intensive,  hence  its  suitability  for  large  scale  simulation  is 
questionable. 

Barbero  and  De  Vivo  (2001)  developed  a  2-D  (plane  stress)  model  for  progressive 
damage  based  on  the  use  of  a  symmetric  second  order  damage  tensor  whose  eigenvalues 
represent  the  density  of  distributed  microcracks  oriented  normal  to  the  principal  material 
coordinate  directions.  This  assumption  regarding  the  form  of  the  damage  tensor  permits 
considerable  simplification  in  the  overall  implementation  of  damage  mechanics  while 
still  maintaining  a  sufficient  level  of  generality  to  adequately  describe  the  observed 
behavior  of  fiber  reinforced  composite  laminates.  In  addition,  damage  evolution  and 
stiffness  reduction  are  computed  for  the  pre- homogenized  composite  material  which 
further  simplifies  the  formulation.  The  model  was  extended  by  Barbero  and  Lonetti 
(2002)  to  include  plasticity,  and  further  extended  by  Lonetti  et  al.  (2003)  to  include 
triaxial  orthotropic  damage  in  terms  of  three  damage  eigenvalues.  Lonetti  et  al.  (2003) 
used  the  triaxial  damage  model  in  conjunction  with  3-D  finite  elements  that  each  span  the 
entire  thickness  of  the  laminate  thus  encompassing  all  of  the  material  layers  in  a  single 
element. 

Robbins  et  al.  (2005)  discussed  the  numerical  treatment  of  3-D  continuum  damage 
mechanics  within  the  context  of  a  macroscopic  finite  element  model  based  on  the  first 
order  shear  deformation  theory  of  laminates.  The  authors  studied  the  effect  of  various 
modeling  parameters  on  predicted  damage  evolution  and  global  failure.  These 
parameters  included  2-D  mesh  density,  2-D  element  type  (i.e.  polynomial  order),  element 
integration  scheme  and  element  distortion  level.  The  progressive  damage  solution  was 
shown  to  exhibit  convergence  with  increasing  mesh  density  and  was  furthermore  shown 
to  exhibit  minimal  sensitivity  to  changes  in  the  remaining  modeling  parameters.  Robbins 
and  Reddy  (2006)  compared  the  damage  predictive  capabilities  of  ESL  models  and 
layerwise  models,  and  concluded  that  ESL  models  consistently  under-predict  local 
damage  evolution  and  consistently  over-predict  global  failure  loads  due  to  the  inability  of 
ESL  models  to  correctly  resolve  the  transverse  shear  stress  maxima. 

2.  Objectives  and  Methods 

While  a  review  of  the  literature  reveals  numerous  studies  devoted  to  finite  element 
modeling  of  the  response  of  actuated  plates,  the  issue  of  structural  integrity  of  such 
composite  systems  has  received  very  little  attention  to  date,  despite  the  fact  that  this 
particular  concern  has  been  raised  in  numerous  papers.  The  material  mismatches  that  are 
inherent  in  actuated  plates  give  rise  to  severe  stress  concentrations  that  are  likely  to  cause 
material  damage,  particularly  in  the  adhesive  layer  used  to  bond  the  actuator  or  sensor  to 
the  structural  substrate.  While  there  have  been  a  few  studies  that  have  investigated  the 
effect  of  damaged  adhesive  bonds  on  the  overall  functionality  of  actuated  plates,  the 
adhesive  damage  itself  was  assumed,  not  predicted  (Sun  et  al.,  2001;  Kim  and  Jones, 
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1996).  To  date,  there  have  been  no  reported  efforts  to  predict  the  mode,  severity  and 
extent  of  adhesive  damage  that  occurs  during  the  application  of  anticipated  service  loads. 

In  response  to  this  deficiency,  the  overall  objective  of  the  present  study  is  to  develop 
and  demonstrate  a  finite  element  methodology  that  is  capable  of  modeling  progressive 
damage  in  adhesive  bond  layers  of  actuated  plates.  This  finite  element  methodology  is 
used  to  investigate  the  mode,  severity  and  extent  of  distributed  microscopic  damage  that 
accumulates  in  the  adhesive  bond  layers  of  a  simply- supported  actuated  plate  that  is 
subjected  to  bending  actuation  and/or  distributed  transverse  loading.  In  addition,  the 
study  quantifies  the  relative  contributions  to  adhesive  damage  from  transverse  shear 
deformation  and  transverse  normal  deformation,  and  further  quantifies  the  effect  of  the 
‘predicted’  adhesive  damage  on  the  overall  functionality  of  the  actuated  plate. 

The  computational  model  that  is  used  to  accomplish  these  goals  is  developed  by 
incorporating  a  3-D  continuum  damage  mechanics  model  into  a  hierarchical, 
displacement -based  finite  element  model  that  is  specifically  formulated  to  facilitate  the 
simulation  of  composite  laminate  behavior.  The  hierarchical  finite  element  formulation 
is  ideally  suited  for  the  present  study  since  it  permits  the  assumed  kinematics  of  the  entire 
model  (or  any  given  element)  to  be  easily  changed.  This  feature  is  particularly  relevant 
for  the  present  focus  problem  where  the  effect  of  kinematic  assumptions  on  predicted 
adhesive  damage  is  not  well  understood.  The  hierarchical  model  includes  the  first-order 
shear  deformation  model  (FSD),  type-I  layerwise  models  (LW1)  and  type- II  layerwise 
models  (LW2)  as  special  cases.  Both  the  LW1  layerwise  model  and  the  more  familiar 
FSD  model  use  a  reduced  constitutive  matrix  that  is  based  on  the  assumption  of  zero 
transverse  normal  stress;  however,  the  LW1  model  also  includes  discrete  layer  transverse 
shear  effects  via  inplane  displacement  components  that  are  C°  continuous  with  respect  to 
the  thickness  coordinate.  The  LW2  layerwise  model  utilizes  a  full  3-D  constitutive 
matrix  and  includes  both  discrete  layer  transverse  shear  effects  and  discrete  layer 
transverse  normal  effects  by  expanding  all  three  displacement  components  as  C° 
continuous  functions  of  the  thickness  coordinate. 

Since  the  prediction  of  local  damage  evolution  requires  the  use  of  a  macroscopic 
model  that  is  capable  of  delivering  accurate  3-D  stress  and  strain  distributions  within 
each  constituent  material  of  the  actuated  plate,  the  use  of  a  first  order  shear  deformation 
(FSD)  model  is  not  appropriate  for  the  task  (Robbins  and  Reddy,  2006);  consequently,  all 
simulations  in  the  present  study  are  performed  using  layerwise  finite  elements. 

The  particular  3-D  continuum  damage  model  used  in  this  study  is  not  intended  to  be 
novel;  rather,  it  is  simply  intended  to  be  representative  of  the  many  published  damage 
models  that  have  been  discussed  in  the  literature  for  predicting  damage  evolution  and 
stiffness  reduction  in  composite  laminates.  For  completeness  sake,  the  development  of 
the  3-D  damage  model  is  fully  described  in  this  paper  and  shares  many  similarities  with 
the  developments  found  in  Barbero  and  De  Vivo  (2001)  and  Lee  et  al.  (1985).  The 
numerical  implementation  of  this  particular  damage  model  is  described  in  detail  by 
Robbins  et  al.  (2005)  and  Robbins  and  Reddy  (2006). 

It  should  be  emphasized  that  the  objective  of  this  problem  is  not  necessarily  to 
produce  a  solution  that  agrees  very  closely  with  experimental  results,  but  rather  to 
demonstrate  that  the  present  methodology  is  adequate  for  such  simulations  and  can  be 
used  to  investigate  the  fundamental  characteristics  of  localized  adhesive  bond  damage 
and  its  impact  on  the  overall  functionality  of  the  actuated  plate.  In  order  to  differentiate 
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the  relative  contributions  to  local  adhesive  damage  from  transverse  shear  deformation 
and  transverse  normal  deformation,  both  the  LW1  and  LW2  layerwise  models  are  used  to 
obtain  solutions.  The  scope  of  the  present  study  is  restricted  to  the  prediction  of 
progressive  damage  and  local  material  failure  that  occur  in  the  adhesive  bond  layers  of 
actuated  plates  under  quasi- static,  monotonic  increasing  load  conditions. 

The  remainder  of  the  paper  is  organized  as  follows.  In  Section  3,  the  3-D  continuum 
damage  mechanics  model  is  developed  for  the  special  case  of  orthotropic  damage, 
culminating  in  the  damaged  constitutive  relations  and  the  governing  equations  that  drive 
the  evolution  of  the  internal  damage  variable.  Section  4  describes  the  hierarchical 
laminate  model  that  is  used  to  produce  the  series  of  macroscopic  models  that  are  used  in 
the  study.  Section  4  also  contains  a  brief  description  of  the  numerical  implementation  of 
the  3-D  continuum  damage  mechanics  model  into  the  hierarchical  macroscopic  finite 
element  model.  In  Section  5,  numerical  examples  are  presented  to  demonstrate  the 
fundamental  characteristics  of  adhesive  damage  predicted  by  the  model.  Finally, 
concluding  remarks  are  given  in  Section  6. 

3.  3-D  continuum  damage  mechanics  model 

On  the  microscopic  scale,  damage  is  characterized  by  molecular  bonds  that  have  been 
broken  and  are  no  longer  available  to  support  loads.  As  damage  accumulates,  the 
molecular  bond  density  of  the  material  decreases,  resulting  in  a  reduction  of  stiffness  and 
consequently  a  nonlinear  (softening)  stress/strain  relationship.  Typically  these  broken 
bonds  are  not  uniformly  distributed  in  the  material,  tending  instead  to  nucleate  into 
microcracks  and  microvoids.  Figure  1  illustrates  the  phenomenological  differences 
between  the  nonlinear  softening  material  behavior  that  occurs  in  the  case  of  progressive 
damage  and  the  more  familiar  case  of  plasticity.  The  plastic  behavior  shown  in  Figure  la 
is  associated  with  the  development  of  slip  planes  that  cause  a  permanent  plastic  strain  ep. 
However,  the  plastic  flow  does  not  affect  the  bond  density  of  the  material,  and 
consequently  the  material  does  not  experience  a  permanent  reduction  in  stiffness.  In 
contrast,  the  nonlinear  behavior  in  Figure  lb  is  due  solely  to  damage  (without  plastic 
flow)  and  thus  exhibits  a  permanent  reduction  in  stiffness.  However,  upon  removal  of  the 
loads,  the  damaged  material  does  not  exhibit  a  permanent  deformation.  Figure  lc  shows 
nonlinear  material  behavior  that  is  caused  by  both  plasticity  and  damage,  and 
consequently  exhibits  both  a  permanent  plastic  strain  and  a  permanent  reduction  in 
stiffness. 
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Figure  1.  Idealized  nonlinear  material  behavior,  (a)  Nonlinearity  caused  by  plasticity, 

(b)  Nonlinearity  caused  by  damage,  and  (c)  Nonlinearity  caused  by  plasticity  and 
damage. 

3.1  Damage  Tensor 

Within  the  context  of  continuum  damage  mechanics,  distributed  microscopic  damage 
can  be  quantified  by  the  use  of  an  appropriate  tensor  field  that  describes  the  orientation 
and  density  of  microcracks  in  the  material.  This  damage  tensor  serves  as  an  evolving 
internal  variable  within  the  framework  of  irreversible  thermodynamics,  and  the  local 
value  of  the  damage  tensor  is  used  to  define  the  damaged  (i.e.  reduced)  stiffness  of  the 
material  Various  forms  of  the  damage  tensor  that  have  been  proposed  in  the  literature 
include  scalars,  vectors,  2nd  order  tensors,  and  4th  order  tensors  (Skrzypek  and 
Ganczarski,  1999;  Lemaitre,  1996).  For  the  sake  of  computational  efficiency,  one  should 
choose  the  simplest  tensorial  form  that  is  capable  of  accurately  describing  the 
distribution,  density,  and  orientation  of  microcracks  and  microvoids  for  the  intended  type 
of  problem.  For  problems  that  can  be  idealized  as  beams,  plates  or  shells,  distributed 
damage  can  be  described  to  an  acceptable  level  of  accuracy  by  a  symmetric  2nd  order 
tensor  D  whose  principal  directions  are  assumed  to  coincide  with  the  principal  material 
directions,  i.e.  orthotropic  damage  (Barbero  and  DeVivo,  2001).  In  the  event  that  a 
particular  material  happens  to  be  isotropic  in  its  undamaged  state,  the  principal  material 
directions  are  simply  chosen  to  coincide  with  the  global  coordinate  system  directions,  i.e., 
two  orthogonal  in-plane  directions  and  the  transverse  direction.  The  eigenvalues  of  the 
orthotropic  damage  tensor  D  (denoted  Di,  D2,  and  D3)  have  a  simple  physical 
interpretation,  namely,  the  ith  eigenvalue  D,  represents  the  effective  fractional  reduction  in 
load  carrying  area  on  planes  that  are  perpendicular  to  the  ith  principal  material  direction 
These  damage  eigenvalues  are  illustrated  in  Figure  2  for  the  special  case  of  a  fiber 
reinforced  composite  material  where  the  first  principal  material  direction  is  chosen  to 
coincide  with  the  fiber  direction. 

The  eigenvalues  of  the  damage  tensor  must  be  in  the  range  0  <  D;  <  1  where  f)  =  0 
corresponds  to  a  complete  lack  of  microcracks  that  are  normal  to  the  ith  principal  material 
direction,  while  Q  =  1  corresponds  to  a  complete  separation  of  the  material  across  planes 
that  are  normal  to  the  ith  principal  material  direction. 
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A)  Distribution  of  B)  Distribution  of  C)  Distribution  of 

microcracks  normal  microcracks  normal  microcracks  normal 

to  the  ‘1’  direction  to  the  ‘2’  direction  to  the  ‘3’  direction 

Figure  2.  Representative  volume  elements  showing  distribution  of  microcracks 
described  by  damage  eigenvalues  Dj,  D2  and  D3.  D,  is  interpreted  as  the  effective 
fractional  reduction  in  load  carrying  area  due  to  microcracks  that  are  normal  to  the  ith 
princ  ipal  material  direction. 

3.2  Damaged  Constitutive  Relations  (Stiffness  Reduction  Scheme) 


In  order  to  relate  the  damage  eigenvalues  to  material  stiffness  reduction,  we  first 
define  the  concept  of  effective  stress,  denoted  <3 .  which  represents  an  attempt  to  assign 
the  actual  internal  material  forces  to  cross  sectional  areas  that  have  been  reduced  due  to 
the  presence  of  damage.  Qualitatively  speaking,  we  expect  the  effective  stress  CJ  to  be 
greater  than  the  apparent  stress  G  (which  uses  the  undamaged  cross  sectional  area); 
however,  the  precise  manner  of  defining  the  effective  stress  is  somewhat  arbitrary,  and 
consequently  numerous  methods  for  defining  effective  stress  have  been  proposed  in  the 
literature  (Skrzypek  and  Ganczarski,  1999).  In  the  present  discussion,  we  use  the 
symmetric  effective  stress  tensor  <3  described  by  Cordebois  and  Sidoroff  (1982)  which 
can  be  expressed  as 


O  =  ^/I-d)  CJ-f-v/1- d)TT  j  =  M_1:C  (1) 

where  M  is  a  4th  order  symmetric  tensor  known  as  the  damage  effects  tensor.  When 
expressed  in  the  principal  material  coordinate  system  using  contracted  notation,  M  can  be 
represented  by  the  following  diagonal  6x6  matrix. 
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In  this  case,  the  effective  stress  components  of  Eq.  (1)  can  now  be  expressed  in  the 
principal  material  coordinate  system  using  contracted  notation  as  follows. 

_  o 

°'  =  (l — d")’  f°r  '=1  -2,3  (no  summation  implied)  (la) 


Vl-D2A/l-£>3  ’  5  V1-  A-nA-A  ’  6 


(lb,c,d) 


Having  adopted  a  definition  for  the  effective  stress  tensor,  we  can  now  define  the 
components  of  the  effective  strain  tensor  8  by  insisting  that  the  same  elastic  strain  energy 
density  should  be  obtained  regardless  whether  it  is  computed  with  the  apparent  stress  <3 
and  the  apparent  strain  8,  or  the  effective  stress  <5  and  effective  strain  8  .  This  is  the  so- 
called  principle  of  equivalent  elastic  strain  energy  density,  which  can  be  expressed  as 
(1/2) a  :e  =  (l/2)a  :e  .  Substitution  of  M :  a  =  a  and  noting  that  a ,  e  and  M  are 
symmetric,  yields  8  =  M:8=8:M.  In  this  case,  the  effective  strain  components  can  be 
expressed  in  the  principal  material  coordinate  system  using  contracted  notation  as 
follows. 


8,  =  e,(l  -  D, ) ,  for  i=  1 ,2,3  (no  summation  implied)  (3a) 

e4  =e4Vl -D2Jl-D3  ,  e5  =85^/l-D3A/l-A  ,  e6  =e6^/l-D1%/l- D2  (3b,c,d) 

The  presence  of  distributed  microscopic  damage,  as  defined  by  the  internal  variable 
D,  affects  the  Helmholtz  free  energy  density  'F,  which  is  postulated  to  be  the  sum  of  an 
elastic  strain  energy  density  cp  and  a  free  energy  T1  ( f3 )  that  is  associated  with  damage 
hardening. 


?  =  ?  (8,  D,  p )  =  — cp(8,D)  +  ?  (p)  (4) 

P 

In  Eq.  (4),  the  internal  variable  (3  is  a  non-dimensional  damage  hardening  parameter  that 
characterizes  the  overall  state  of  damage  and  is  used  to  control  the  process  of  damage 
hardening.  Assuming  that  the  homogenized  composite  material  exhibits  a  linear  elastic 
stress/strain  relationship  up  to  the  point  of  damage  initiation  or  up  to  the  point  of  further 
damage  progression,  the  elastic  strain  energy  density  can  be  expressed  as 
cp  (8.  D)  =  (1/2)0 : 8  =  (1/2)8  :  C  :  8 ,  where  C  =  C(D)  is  the  4th  order  damaged  material 

elasticity  tensor.  The  dependence  of  C  on  the  damage  tensor  D  is  defined  by  requiring 
that  the  elastic  strain  energy  density  of  the  damaged  material  that  is  subjected  to  the 
apparent  strain  8  should  be  equivalent  to  the  stain  energy  density  of  the  undamaged 
material  that  is  subjected  to  the  effective  (or  reduced)  strain  8  .  Upon  introducing  the 
expression  (8  =  M:8=8:M)  into  the  strain  energy  equivalence  principle,  the  damaged 

material  elasticity  tensor  can  be  defined  as  C  =  M:C:M.  C  can  be  more  simply 
expressed  in  the  principal  material  coordinate  system  using  contracted  notation  as  a 
symmetric  6x6  matrix  that  is  indicative  of  an  orthotropically  damaged  material.  In  this 
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case,  the  nonzero  components  of  C  are  given  by  Eq.  (5)  where  repeated  indices  do  not 
imply  summation. 


Op  =Cap(l-DJ(l-Dp)  fora, (3=1, 2, 3  (5a) 

C44  =  C44(1-D2)(1-D3),  C55  =  c55(l- )(!-£>,),  C66  =C66(1-£>1)(1-D2)  (5b,c,d) 


3.3  77ze  damage  surface  and  the  damage  evolution  equations 


After  having  defined  stiffness  reduction  in  terms  of  the  damage  eigenvalues,  we  must 
now  develop  the  evolution  equations  for  the  damage  eigenvalues.  Increments  in  the  free 
energy  density  p*P  can  be  expressed  in  terms  of  increments  in  the  kinematic  variables  and 
internal  variables  (£,  D,  (3)  as 


A(p?)  = 


*Pl2:Ae  +*£!2:ad 

38  3D  3(3 


(6) 


where  the  partial  derivatives  in  Eq.  (6)  define  the  generalized  thermodynamic  forces  that 
are  energy  conjugate  to  the  variables  £,  D  and  (3.  Partial  differentiation  of  p'F  with 
respect  to  the  apparent  strain  tensor  £  defines  the  apparent  Cauchy  stress  tensor  o  that 
acts  on  the  damaged  material. 

?( 1  “  ^ 

3(p?)  _  3cp 


o  = 


— £:C:£ 
2 


3e 


3e 


3e 


=  C:£ 


(7) 


Partial  differentiation  of  p'P  with  respect  to  the  damage  tensor  D  yields  a  2  order, 
symmetric,  negative  semidefinite  tensor  Y  that  is  energy  conjugate  to  the  damage  tensor. 
However,  it  is  common  convention  to  include  a  sign  change  in  the  definition,  i.e.  Y  = 

— 3  ( p  )/3 D .  in  which  case  Y  is  interpreted  as  the  damage  energy  release  rate  tensor. 


3(p?) 

3D 


3cp 

3D 


1 

— £:C:£ 


3D 


3 


1  3C 

-£: - : 

2  3D 


(8) 


Eq.  (8)  can  be  conveniently  expressed  in  the  principal  material  coordinate  system  via  use 
of  Eqs.  (5),  yielding  the  following  expressions  for  the  damage  energy  release  rate 
eigenvalues. 


Y]  —  Cj  j  (1  —  D,  )£  j  +  C12  (1—  D2  )£  j£t  +  C13(l —  D 3  )£  j£  3  +  ■ 

^2  —  C21  (1  -Di )£ j£ 2  +  c22(l  —  D2)e2  C23 (1  —  ^-^3 )£ 2^3 

F3  =  C3 j (1 —  Z)j )£  j£ 3  +  (1 —  Z)9  )£  3£  9  +  C33 (1 —  Z)3 )£ 3  + 


^(1-D,)£52+^(1-D2)£62 

^-(1-D3)£42^(1-A)£62 

^i(l-D2)£42+-^(l-D,)£2 


(9a,b,c) 
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Partial  differentiation  of  p'P  with  respect  to  the  overall  damage  parameter  (3  yields  a 
generalized  thermodynamic  force  referred  to  as  the  damage  hardening  variable  y(|3). 


Y  =Y(P)  =  ^(£72  =  an 
ap  ap 

A  linear  relationship  y  =  c i|3  is  often  used  for  damage  hardening  (Voyiadjis  and  Kattan, 
1999;  Lee  et  aL,  1985);  however,  an  exponential  form  such  as  y  =  q( \-e^^2 )  is  often 

better  able  to  represent  the  damage  hardening  characteristics  observed  in  fiber- reinforced 
composites  (Barbero  and  DeVivo,  2001). 

It  is  postulated  that  damaging  behavior  can  be  distinguished  from  non-damaging 
behavior  on  a  local  basis  by  a  convex  damage  surface  g(  Y,y)  =  0  that  is  defined  in  energy 
release  rate  space,  where  g(Y,y)<0  indicates  a  non-damaging  state,  g(Y,y)=0  indicates  a 
damage  inducing  state,  and  g(Y,y)>0  is  understood  to  be  inadmissible.  The  damage 

1  8  90  99 

surface  is  often  assumed  to  be  a  quadratic  function  of  the  energy  release  rate  tensor  ’  ’ 
As  a  specific  example,  we  consider  the  following  damage  surface,  expressed  in  the 
principal  material  coordinate  system. 


g( Y,y)=  V/11f12  +  722f22  +  733f32  -  (y0+y(P)) 


In  Eq.  (11),  Jn,  J22,  and  J33  represent  experimentally  derived  material  constants  that 
define  the  damage  tolerance  of  the  material,  and  effectively  control  the  shape  of  the 
damage  surface.  The  expression  yo  +  y([3)  defines  the  current  damage  threshold  and  thus 
controls  the  size  of  the  damage  surface.  The  material  constant  yo  defines  the  initial 
damage  threshold  of  the  un-damaged  material,  and  y( [3 )  provides  the  increase  in  the 
damage  threshold  (i.e.  damage  hardening)  that  occurs  as  the  damage  tensor  D  evolves. 
Note  that  the  energy  release  rate  eigenvalues,  as  defined  in  Eqs.  (9),  are  positive 
semidefinite  quantities,  regardless  whether  the  material  exhibits  tensile  deformation  or 
compressive  deformation.  Thus  the  simple  damage  surface  defined  by  Eq.  (11)  is  unable 
to  account  for  materials  that  are  more  easily  damaged  in  tension  than  compression.  For  a 
discussion  of  models  that  incorporate  this  unilateral  effect,  the  reader  is  referred  to 
Ladeveze  (2002)  and  Dragon  (2002). 

The  damage  surface  of  Eq.  (11)  defines  the  limits  of  non-damaging  behavior  at  an 
arbitrary  point  in  the  material.  As  long  as  the  material  exhibits  variables  (Yi,  Y2,  Y3,  y) 
that  identify  points  interior  to  the  current  damage  surface  g(Y,y)  =  0,  then  no  further 
damage  occurs  at  the  local  material  point  where  Yi,  Y2,  Y3,  and  y  are  measured. 

However,  when  these  variables  cause  the  damage  surface  g(  Y,y)  to  be  reached  or 
exceeded,  then  damage  is  assumed  to  progress  in  a  manner  that  is  consistent  with  the 
Kuhn-Tucker  conditions  (g  <  0,  dy  >  0,  and  g-dy  =  0)  which  imply  that  the  condition  g=0 
(or  alternately  dg  =  0)  must  be  maintained  during  a  damaging  process.  In  other  words, 
any  attempt  to  reach  a  point  (Y,y)  that  would  cause  g(  Y,y)>0  results  in  a  simultaneous 
increase  in  the  damage  tensor  D  and  the  damage  hardening  variable  y  in  such  a  way  as  to 
maintain  the  condition  g(Y,y)=0  throughout  the  damaging  process. 
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The  exact  manner  in  which  the  internal  variables  D  and  y  evolve  is  provided  by  the 
Principle  of  Maximum  Dissipation  which  states  that  during  a  damaging  (i.e.  dissipative) 
process,  the  actual  values  of  the  thermodynamic  forces  Y  and  y  will  cause  the  dissipation 
power  density  $  to  be  extremized,  subject  to  the  constraint  that  g=0  (Lee  et  al.,  1985). 
For  the  case  of  elastic  damage  (i.e.  damage  without  accompanying  plasticity),  the 
dissipation  power  density  can  be  expressed  in  the  principal  material  coordinate  system  as 

<F(Y,y)  =  F/D/-yp>0  (12) 

In  order  to  extremize  the  dissipation  power  density  subject  to  the  constraint  g(Y,y)=0,  we 
introduce  the  Lagrange  multiplier  X  and  form  the  modified  objective  function  <f>  . 


<D*(Y,y)  =  ¥,£>,-  yp  -Xg  (13) 

Extremization  of  <f>  with  respect  to  the  variables  Y  and  y  yields 


dd> 

dy 


Dj  -  X-^-  =  0 

dr, 

-p  -X^-  =0 

dy 


-»  bj=x^- 

dY, 

->  P  =-^ 

dy 


(14a) 

(14b) 


Eqs.  (14)  express  the  rate  of  change  of  the  damage  tensor  components  and  the  overall 
damage  parameter  in  terms  of  the  common  Lagrange  multiplier  X.  The  value  of  X  can  be 
determined  from  the  consistency  condition  which  requires  that  during  a  damaging 
process,  the  damage  surface  should  evolve  in  such  a  way  that  g=0  (or  alternately  g  -  0 )  is 
maintained. 


g(  Y,y) 


_  dg 

d  r, 


dg 

Y,  +yr~i  ~  0 

dy 


(15) 


Eqs.  (14)  and  (15)  provide  the  relationships  that  are  needed  to  determine  the  evolution  of 
the  damage  tensor  D  and  the  overall  damage  parameter  (1 .  The  Lagrange  multiplier  ?  can 
be  eliminated  from  this  system  of  equations  to  yield 


h  ) 

f  dg  dg  } 

d  Y 

U  1  K 

^  _ 

dr,  dr. 

dg  dy  dg 

dr, 

dg  dry  dg 

v  dy  dp  dy  y 

v  dy  5P  dy  y 

(16a) 
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3s  t' ) 

(  as 

dYK 

dg  _  dYK 

dg  dy  dg 

dy  dg  dy 

dy  dp  dy 

1  dy  ap 

WkYk 


(16b) 


where  the  quantities  Z [k  and  Wk  have  been  introduced  for  convenience.  The  time 
variable  can  be  eliminated  from  Eq.  (16)  to  yield  an  expression  which  relates  increments 
in  the  damage  eigenvalues  and  overall  damage  parameter  to  increments  in  the  energy 
release  rate  eigenvalue. 


dDi  =  Zik  dYic, 


WK  dYK 


(17a,b) 


In  its  present  form,  Eq.  (17)  is  not  particularly  useful  since  it  is  not  practical  to  impose 
increments  in  the  energy  release  rate  eigenvalues.  However,  Eq.  (17)  can  be  converted  to 
an  imposed  strain  form  by  invoking  Eq.  (9)  [Y  =  Y(e  ,D)]  to  express  dYK  as 


?)Y,  })Yk 

d  Y K  —  — dD  n  +  — — da 

K  dDQ  Q  <ka  “ 


where  we  have  adopted  contracted  notation  for  the  strain  tensor  components  (i.e.  £i  =£n, 
£2  =£22,  £3  =£33,  £4  =  2£23,  £5  =  2£3i,  £6  =  2e  12).  Substitution  of  Eq.  (18)  into  Eq.  (17) 
permits  the  governing  differential  equations  of  damage  evolution  to  be  written  in  terms  of 
dD,  dp  and  d£. 


dYK  dYK 

dD j  -  Z1K  ——dD  +  ZIK  ——da( 
dDQ  y  dea 

dV=WK-^dDQ  +  WKj^d£a 
dD  y  dea 


(19a) 

(19b) 


Equation  (19)  can  be  solved  for  d Di  (1=1, 2, 3)  and  c/[).  yielding  the  following  imposed 
strain  form  of  the  damage  evolution  equations. 


=  [TKrfe}. 


In  Eq.  (20),  the  4x6  matrix  [T]  can  be  expressed  as  [A]"1  [B],  where  the  4x4  matrix  [A] 
and  the  4x6  matrix  [B]  are  defined  as 
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[A]  = 


1-Z 


IK 


-z 


2K 


-z 


2K 


-Ww 


ay, 

:  BD, 

„  37a. 

7  K 

“  sa 

7  3Yk 

IK  BD, 

0 
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a  d, 

l-Z  ^ 
“3D, 

7  3YK 

2K  dD 3 

0 

ay, 

BD, 

7  3Y, 

M  BO, 

1  Z3  KdYR 

3K  dD 3 

0 

ay, 

3D, 

3  D2 

-vy,  ^ 

A  3D, 

1 

[B]  = 


Z 

Z 

z 


K 


'1Jr  38, 


a: 


3^ 

'2*  3e, 

3F, 

3*  3e, 


K 


3  Y, 


K 


3e, 


7 

^ia: 
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< 
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an 
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an 
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3n 
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7 
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7 
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7 
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ay. 

WK 

ay. 
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(21a) 


(21b) 


5.4  Example  -  Damaged  Response  of  the  Adhesive  Material 

To  demonstrate  the  material  response  characteristics  predicted  by  the  3-D  continuum 
damage  mechanics  model  we  consider  the  homogeneous  deformation  and  loading  of  a 
brittle  epoxy  adhesive  material  that  is  typical  of  the  adhesives  used  to  bond  actuators  to 
structural  substrates.  In  its  un-damaged  (virgin)  state,  the  adhesive  material  is  assumed  to 
be  isotopic,  with  modulus  E  =  6.8  GPa  and  Poisson’s  ratio  v  =  0.3.  The  damage  surface 
constants  and  the  damage  hardening  constants  for  the  adhesive  material  are  assumed  as 
Ji  i  =  J22  =  J33  =  1,  Y(P)  =  2.5(  106)p.  and  y(0)  =  0.  Two  different  load  cases  are 
considered.  The  first  load  case  involves  imposing  a  homogeneous  state  of  simple  shear 
deformation  in  the  1-3  plane  (i.e.,  85 >0.  all  other  e«  are  zero).  This  particular  mode  of 
deformation  contributes  equally  to  the  evolution  of  damage  eigenvalues  Di  and  D3,  thus 
resulting  in  the  accumulation  of  microcracks  that  are  oriented  normal  to  the  T’  and  ‘3’ 
directions  respectively.  Part  A  of  Figure  3  shows  the  non-linear  softening  stress/strain 
relationship  predicted  by  the  damage  model.  Note  that  the  maximum  value  of  shear 
stress  attained  in  this  displacement -controlled  test  is  approximately  36  MPa.  For  a 
similar  load- controlled  shear  test,  the  adhesive  material  simply  exhibits  complete  material 
failure  during  any  attempt  to  increase  the  imposed  shear  stress  beyond  36  MPa.  Part  B  of 
Figure  3  shows  the  evolution  of  damage  eigenvalues  Di  and  D3  caused  by  the  imposed 
shear  deformation.  Examination  of  Parts  A  and  B  reveals  that  when  the  maximum  shear 
stress  is  attained,  the  magnitude  of  the  damage  eigenvalues  Di  and  D3  is  approximately 
0.24. 

Next,  consider  the  case  where  the  adhesive  material  is  subjected  to  uniaxial  extension 
in  the  ‘3’  direction.  During  the  imposed  extension,  the  material  is  free  to  undergo  lateral 
contraction  in  the  1  and  2  directions  so  that  it  maintains  a  uniaxial  stress  state.  This 
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particular  mode  of  deformation  contributes  only  to  damage  eigenvalue  D3,  resulting  in 
the  accumulation  of  microcracks  that  are  oriented  normal  to  the  ‘3’  direction.  Part  C  of 
Figure  3  shows  the  damaged  uniaxial  stress/strain  relationship  predicted  by  the  damage 
model.  Note  that  the  maximum  value  of  uniaxial  normal  stress  that  can  be  supported  by 
the  adhesive  material  is  approximately  41  MPa;  any  attempt  to  load  the  adhesive  material 
beyond  this  value  results  in  complete  failure  of  the  material.  Part  D  of  Figure  3  shows 
the  evolution  of  damage  eigenvalue  Eh  caused  by  the  imposed  axial  extension. 
Examination  of  Parts  C  and  D  reveals  that  when  the  maximum  uniaxial  stress  is  attained, 
the  magnitude  of  the  damage  eigenvalue  D3  is  approximately  0.24. 


A)  imposed  shear  deformation 


C)  imposed  uniaxial  extension 


0  0.005  0.01  0.015  0.02  0.025  0.03 

shear  strain 


B)  imposed  shear  deformation 


0  0.005  0.01  0.015 

uniaxial  normal  strain 
D)  imposed  uniaxial  extension 


Figure  3.  Progressive  damage  in  an  adhesive  material  under  homogeneous  deformation 
and  loading.  A)  Damaged  stress/strain  curve  during  imposed  simple  shear  deformation 
in  the  1-3  plane,  B)  Damage  evolution  during  imposed  simple  shear  deformation  in  the 
1-3  plane,  C)  Damaged  stress/strain  curve  during  imposed  uniaxial  extension  in  the  3 
direction,  D)  Damage  evolution  during  imposed  uniaxial  extension  in  the  3  direction. 


In  Section  5,  this  particular  adhesive  material  description  is  used  in  modeling  the 
progressive  adhesive  damage  that  occurs  in  a  simply- supported  actuated  plate  subjected 
to  bending  actuation  and  distributed  transverse  loading. 
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4.  Hierarchical,  variable -kinematic,  finite  element  model 

To  facilitate  the  goals  of  the  present  study,  a  2-D,  hierarchical,  displacement -based, 
variable -kinematic  finite  element  model  is  used  for  all  simulations  that  are  performed.  The 
main  attribute  of  the  variable -kinematic  finite  element  model  that  makes  it  suitable  for  the 
present  study  is  that  the  kinematics  and  constitutive  relations  of  any  given  element  can  be 
conveniently  changed,  thus  allowing  any  given  element  to  represent  a  wide  variety  of 
different  laminate  models  ranging  from  the  very  simple  to  the  very  complex.  In  addition, 
the  variable-kinematic  finite  element  model  permits  different  types  of  laminate  elements 
(representing  different  types  of  laminate  models)  to  be  conveniently  connected  together  in 
the  same  computational  domain,  thus  permitting  different  subregions  to  be  described  by 
different  mathematical  models. 

Within  the  context  of  the  present  study,  the  2-D,  hierarchical,  displacement-based, 
variable -kinematic  finite  element  is  developed  by  expressing  the  total  displacement  field  as 
the  sum  of  a  low  order  primary  displacement  field  and  a  higher  order  secondary 
displacement  field.  The  primary  displacement  field  is  present  in  all  variable -kinematic 
elements  at  all  times.  The  individual  terms  of  the  secondary  displacement  field  then  serve  as 
relative  displacements  that  can  be  added  to  the  element’s  primary  field  to  provide  higher 
order  kinematics  as  needed.  The  total  displacement  field  is  expressed  as 

Ua(x,y,z)  =  uaFSD(x,y,z)  +  u«LW(x,y,z),  a  =  1,2,3  (22) 

where  ui,  u2,  and  U3  are  the  total  displacement  components  in  the  x,  y,  and  z  directions 
respectively.  The  primary  displacement  field  is  provided  by  nzFSD(x,y,z)  which  represents 
the  assumed  displacement  field  for  the  first  order  shear  deformation  theory  (FSD)  and  is 
expressed  as 

uiFSD(x,y,z)  =  uo(x,y)  +  20x(x,y),  (23a) 

u2FSD(x,y,z)  =  v0(x,y)  +  z0y(x,y),  (23b) 

u3FSD(x,y,z)  =  w0(x,y).  (23c) 

In  Eqs.  (23),  uo(x,y),  vo(x,y)  and  wo(x,y)  represent  the  displacement  of  points  on  the  plate’s 
reference  surface  which  is  normally  chosen  to  coincide  with  the  plate’s  mid -surface.  The 
terms  0x(x,y)  and  0y(x,y)  represent  the  rotation  of  the  inextensible  transverse  normal  fiber  in 
the  xz  and  yz  planes  respectively.  The  independent  rotations  0x(x,y)  and  0y(x,y)  are  not 
required  to  be  equal  in  magnitude  to  the  slopes  dwo/dx  and  dwo/dy,  thus  the  FSD 
displacement  field  includes  a  rudimentary  transverse  shear  strain  that  is  constant  through  the 
thickness  of  the  laminate.  Since  the  FSD  displacement  field  does  not  explicitly  include 
transverse  normal  strain,  it  is  intended  to  be  used  in  conjunction  with  a  reduced  constitutive 
matrix  that  is  based  on  the  assumption  of  zero  transverse  normal  stress. 

The  secondary  displacement  field  in  Eq.  (22),  labeled  as  UaLW(x,y,z),  represents  the 
assumed  displacement  field  for  a  full  3-D  layerwise  theory28"30  which  is  characterized  by 
displacement  components  that  are  piecewise  continuous  (specifically,  C°  continuous)  with 
respect  to  the  thickness  coordinate.  The  layerwise  displacement  field  is  included  as  an 

ccr\ 

optional,  incremental  enhancement  to  the  primary  displacement  field,  Ua  (x,y,z),  so  that 
the  element  may  have  full  or  partial  3-D  modeling  capability  when  needed.  Depending  on 


16 


the  desired  level  of  accuracy,  the  variable -kinematic  element  can  use  none,  part,  or  all  of  the 
layerwise  field  u«LW(x,y,z)  to  create  a  series  of  different  elements  having  a  wide  range  of 
kinematic  complexity.  For  example,  discrete  layer  transverse  shear  effects  can  be  added  to 
the  element  by  including  uiLW(x,y,z)  and  U2LW(x,y,z).  Discrete  layer  transverse  normal 
effects  can  be  added  to  the  element  by  including  U3LW(x,y,z).  The  layerwise  field  can  be 
expressed  as 

uiLW(x,y,z)  =  Uj(x,y)cpj(z) 
u2LW(x,y,z)  =  Vj(x,y)cpj(z) 
u3LW(x,y,z)  =  Wj(x,y)cpj(z) 

where  the  repeated  subscript  j  implies  summation  over  j=l,2,...,n.  The  functions  cpj(z) 
(]=l,2,..,n)  are  1-D  Lagrangian  interpolation  functions  associated  with  n  nodes  distributed 
through  the  laminate  thickness,  located  at  Zj  (j=l,2,..,n).  Thus  the  through- the-thickness 
variation  of  the  displacement  components  is  defined  in  terms  of  a  1-D  finite  element 
representation  with  C°  continuity  of  the  interpolants.  The  1-D  interpolants  Uj(x,y),  Vj(x,y), 
and  Wj(x,y)  represent  additions  to  the  displacement  components  ui,  U2,  and  u3  on  the  planes 
defined  by  z=Zj  (j=l,2,..,n).  It  should  be  noted  that  the  layerwise  field  given  by  Eqs.  (24)  is 
sufficiently  general  to  model  any  of  the  deformation  modes  that  can  be  modeled  by  the  FSD 
field  given  in  Eqs.  (23);  thus,  for  elements  that  use  all  the  variables  shown  in  Eqs.  (23)  and 
(24),  there  will  be  five  redundant  variables  that  must  be  set  to  zero  (or  ignored)  to  permit  a 
unique  solution  for  the  remaining  variables.  The  presence  of  the  FSD  variables  is  essential 
for  connecting  different  types  of  elements,  thus  five  of  the  layerwise  variables  should  be  set 
to  zero  (for  example,  Ui=Un=0,  Vi=Vn=0,  Wi=0). 

A  hierarchy  of  three  distinctly  different  types  of  laminate  elements  can  be  obtained  from 
the  composite  displacement  field  of  Eq.  (22),  where  the  individual  displacement  expansions 
are  defined  by  Eqs.  (23)  and  (24).  The  first  and  simplest  type  of  element  is  the  first  order 
shear  deformation  element  (or  FSD  element).  This  element  is  formed  using  Eqs.  (23)  while 
ignoring  Eqs.  (24).  The  FSD  element  uses  a  reduced  constitutive  matrix  that  is  derived 
based  on  the  assumption  of  zero  transverse  normal  stress.  The  second  type  of  element  is  the 
Type- 1  layerwise  element  (or  LW1  element).  The  LW1  element  is  formed  using  Eqs.  (23), 
(24a)  and  (24b),  while  ignoring  (24c),  thus  the  LW1  element  includes  discrete  layer 
transverse  shear  effects  but  neglects  transverse  normal  effects  and  consequently  uses  a 
reduced  stiffness  matrix  similar  to  the  FSD  element.  Due  to  the  inclusion  of  discrete  layer 
transverse  shear  effects,  the  LW 1  element  is  applicable  to  thick  laminates  and  often  yields 
results  comparable  to  3-D  finite  elements  while  using  approximately  two  thirds  the  number 
of  degrees  of  freedom.  The  third  and  most  complex  element  is  the  Type-II  layerwise 
element  (or  LW2  element).  The  LW2  element  is  formed  using  both  Eqs.  (23)  and  (24),  thus 
it  is  a  full  3-D  layerwise  element  that  explicitly  accounts  for  all  six  strain  components,  and 
consequently  uses  a  full  3-D  constitutive  matrix.  The  inclusion  of  the  full  layerwise  field 
provides  the  LW2  element  with  both  discrete  layer  transverse  shear  effects  and  discrete 
layer  transverse  normal  effects.  In  terms  of  interpolation  capability  and  number  of  degrees 
of  freedom,  the  2-D  LW2  element  is  equivalent  to  an  entire  stack  of  conventional  3-D  finite 
elements. 

Since  the  first  order  shear  deformation  field  is  present  in  each  of  these  three  different 
types  of  elements,  displacement  continuity  can  be  maintained  throughout  the  mesh,  even  if 


(24a) 

(24b) 

(24c) 
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all  three  element  types  are  simultaneously  present  in  different  parts  of  the  mesh.  This  is 
achieved  by  simply  enforcing  homogeneous  boundary  conditions  on  some  or  all  of  the 
incremental  layerwise  variables  along  the  boundaries  that  separate  incompatible  subregions, 
a  process  that  is  easily  automated.  A  significant  advantage  afforded  by  the  hierarchical 
formulation  is  that  once  the  in-plane  mesh  is  defined,  the  user  can  then  assign  any  of  the 
three  element  types  (FSD,  LW1,  LW2)  to  any  of  the  elements  in  the  2-D  mesh.  Subsequent 
changes  in  the  type  of  any  single  element  or  group  of  elements  can  be  performed  with 
minimal  effort,  thus  facilitating  the  goals  of  the  present  study. 

4.1  Implementation  of  3-D  damage  mechanics  into  the  multilayer  VKFE  element 

The  3-D  continuum  damage  mechanics  formulation  described  in  Section  3  introduces 
two  internal  variables  (the  damage  tensor  D  and  the  overall  damage  parameter  (3)  whose 
evolution  equations  (Eq.  20)  are  developed  within  the  framework  of  irreversible 
thermodynamics.  In  order  to  efficiently  apply  the  damage  formulation  to  fiber- reinforced 
laminated  composite  materials,  several  simplifying  assumptions  were  invoked.  First,  the 
damage  evolutions  are  applied  to  the  homogenized  description  of  the  composite  material 
instead  of  being  separately  applied  to  each  of  its  constituent  materials.  Second,  the 
homogenized  composite  material  is  assumed  to  exhibit  orthotropic  damage,  hence  the 
state  of  damage  can  be  described  by  the  three  eigenvalues  of  a  symmetric  second  order 
damage  tensor  (D  — >  Di,  D2,  D3).  Third,  the  homogenized  composite  material  is  assumed 
to  exhibit  isotropic  damage  hardening  which  can  be  described  by  a  scalar  function  of  the 
overall  damaged  state. 

The  numerical  implementation  of  this  damage  model  in  a  first  order  shear  deformable 
element  was  discussed  by  Robbins  et  al.27  and  is  briefly  summarized  as  follows.  First, 
the  polynomial  order  of  the  2-D  elements  is  restricted  to  linear  or  quadratic  forms.  Finear 
2-D  elements  permit  a  better  representation  of  damage  localization,  while  quadratic 
elements  offer  superior  performance  in  bending  dominated  problems.  Second,  the 
damage  evolution  equations  are  solved  only  at  the  discrete  damage  sampling  points  that 
are  chosen  to  coincide  with  the  reduced  Gaussian  integration  points  of  the  2-D  element. 
Third,  based  on  the  damage  eigenvalues  computed  at  the  damage  sampling  points,  the 
damage  eigenvalues  are  interpolated  over  the  2-D  element’s  reference  surface  using 
polynomials  that  are  one  order  less  than  the  displacement  field.  Fourth,  integration  in  the 
thickness  direction  is  performed  using  full  quadrature ;  three  Gauss  points  per  material 
layer  provides  exact  thickness  integration  of  all  terms  even  if  the  damage  eigenvalues 
exhibit  quadratic  or  cubic  variation  through  the  thickness  of  each  material  layer.  Finally, 
the  2-D  FSD  element  uses  a  reduced  constitutive  matrix  that  is  stored  as  a  full  6x6  matrix 
where  the  transverse  normal  components  (row  3  and  column  3)  are  set  to  zero;  therefore, 
the  2-D  FSD  element  can  directly  utilize  the  full  3-D  damage  mechanics  equations  in 
unaltered  form.  These  implementation  guidelines  are  also  directly  applicable  to 
layerwise  elements  that  use  either  linear  or  quadratic  layers,  and  are  used  in  the  present 
study  for  all  three  element  types  (FSD,  FW1  and  FW2)  that  are  derivable  from  the 
hierarchical  VKFE  finite  element  model. 
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5. 


Numerical  Results 


In  this  section,  the  3-D  progressive  damage  model  is  used  in  conjunction  with  the  2- 
D,  hierarchical,  displacement -based,  variable -kinematic  finite  element  model  to  simulate 
damage  accumulation  in  the  adhesive  bond  layers  of  a  simply- supported,  actuated  plate 
that  is  subjected  to  bending  actuation  and/or  distributed  transverse  loading.  Since  the 
prediction  of  local  damage  evolution  requires  the  use  of  a  macroscopic  model  that  is 
capable  of  delivering  accurate  3-D  stress  and  strain  distributions  within  each  constituent 
material  of  the  actuated  plate,  the  use  of  a  first  order  shear  deformation  (FSD)  model  is 
not  appropriate  for  the  task;  consequently,  all  simulations  are  performed  using  layerwise 
finite  elements.  It  should  be  emphasized  that  the  objective  of  this  problem  is  not 
necessarily  to  produce  a  solution  that  agrees  very  closely  with  experimental  results,  but 
rather  to  demonstrate  that  the  present  methodology  is  adequate  for  such  simulations  and 
to  investigate  the  fundamental  characteristics  of  localized  adhesive  bond  damage  and  its 
impact  on  the  overall  functionality  of  the  actuated  plate.  In  order  to  differentiate  the 
relative  contributions  to  local  adhesive  damage  from  transverse  shear  deformation  and 
transverse  normal  deformation,  both  the  LW1  and  LW2  layerwise  models  are  used  to 
obtain  solutions. 

Figure  4  shows  the  geometry  and  boundary  conditions  of  the  simply- supported 
actuated  plate  that  is  used  throughout  the  present  study.  The  test  specimen  is  composed 
of  a  square  aluminum  plate  and  a  symmetric  pair  of  square  piezoceramic  actuators  that 
are  bonded  to  the  upper  and  lower  surfaces  of  the  aluminum  plate  via  a  thin  adhesive 
layer.  The  relevant  elastic  and  piezoelastic  material  properties  that  are  assumed  for  the 
aluminum  substrate,  the  piezoceramic  actuator  and  the  adhesive  are  listed  below.  Note 
that  the  stiffness  of  the  adhesive  material  is  one-tenth  the  stiffness  of  the  piezoceramic 
actuator. 

Elastic  and  piezoelectric  constants  for  the  un-damaged  isotropic  materials'. 

Aluminum:  E  =  70  GPa,  v  =  0.3, 

Adhesive:  E  =  6.3  GPa,  v  =  0.3, 

Piezoceramic:  E  =  63  GPa,  v  =  0.3,  d3i  =  d32  =  -3.74537(10"7)  mm/volt. 

The  length  of  the  aluminum  plate  (2L)  is  fixed  at  twice  the  length  of  the  piezoceramic 

actuators  (2P).  The  thickness  of  each  piezoceramic  actuator  is  denoted  as  h.  The 
thickness  of  the  aluminum  substrate  is  denoted  as  It,  and  the  thickness  of  the  adhesive 
layer  is  denoted  as  ha.  The  total  thickness  of  the  actuated  region  is  then  H=2/z+2ha+hp. 
While  the  total  thickness  H  of  the  actuated  region  is  varied  in  the  study,  the  ratio  of 
actuator  thickness  to  aluminum  thickness  is  fixed  at  /z/hs=0.25,  and  the  ratio  of  adhesive 
thickness  to  actuator  thickness  is  fixed  at  fi//i=0.10.  The  edges  of  the  aluminum  plate  at 
x=±L  and  y=±L  are  simply  supported.  Due  to  symmetries  that  exist  along  the  x  and  y 
axes,  the  2-D  computational  domain  is  reduced  to  one  quadrant  of  the  actuated  plate. 

Parts  B  and  C  of  Figure  4  show  two  different  thickness  configurations  that  are  considered 
in  the  present  study,  namely,  a  relatively  thick  configuration  characterized  by  an  actuated 
span- to- thickness  ratio  2P/H=2  and  an  aluminum  plate  span-to-thickness  ratio  2L/hp=6, 
and  a  more  moderate  thickness  configuration  characterized  by  an  actuated  span-to- 
thickness  ratio  2P/H=8  and  an  aluminum  plate  span-to-thickness  ratio  2L/hp=24. 
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Two  forms  of  loading  are  considered  for  the  simply- supported  actuated  plate.  The 
first  form  of  loading  consists  of  applying  opposite  electric  fields  to  the  top  and  bottom 
actuators,  thus  inducing  global  bending  of  the  actuated  plate.  The  second  form  of  loading 
consists  of  applying  a  uniform  distributed  transverse  load  to  the  entire  upper  or  lower 
surface  of  the  actuated  plate. 


Simply  Supported  Boundary  Conditions 

At  x=0:  u(0,y,z)=0 

At  x=L:  v(L,y,z)=0,  w(L,y,z)=0 

At  y=0:  v(x,0,z)=0 

At  y=L:  u(x,L,z)=0,  w(x,L,z)=0 

piezoceramic  actuator  thickness  h 
adhesive  bond  thickness  h 

d 

aluminum  plate  thickness  hp 
total  thickness  of  actuated  region  H 

h  =  0.25hp  ha  =  0.10/i 


Figure  4.  A)  View  of  specimen  geometry  in  XY  plane  and  simply  supported  boundary 
conditions,  B)  Thickness  configuration  characterized  by  2P/H=2  and  2L/hp=6,  C) 
Thickness  configuration  characterized  by  2P/H=8  and  2L/hp=24. 

5. 1  Effect  of  mesh  density  and  model  type  on  local  adhesive  layer  deformation 

Since  local  damage  initiation  and  local  damage  progression  are  dependent  on  the 
local  state  of  deformation,  it  is  first  necessary  to  determine  the  effects  of  2-D  mesh 
density  and  mathematical  model  type  on  the  predicted  local  deformation  in  the  adhesive 
bond  layer.  To  this  end,  the  simply- supported  actuated  plate  is  subjected  to  bending 
actuation  without  any  distributed  transverse  loading,  and  linear  elastic  (non  damaging) 
solutions  are  computed  using  the  LW1  and  LW2  layerwise  models  in  conjunction  with 
five  different  levels  of  2-D  mesh  density.  All  five  2-D  meshes  are  composed  of  8-node, 
quadratic  quadrilateral  elements.  Three  of  the  2-D  meshes  are  uniform,  while  the 
remaining  two  meshes  are  nonuniform  (i.e.  graded)  where  the  element  size  decreases  as 
the  two  free  edges  of  the  actuator  are  approached.  The  three  uniform  2-D  meshes 
include  a  6x6  uniform  mesh,  a  12x12  uniform  mesh,  and  a  16x16  uniform  mesh,  denoted 
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respectively  as  6x6u,  12xl2u  and  16xl6u  Since  the  actuated  region  occupies  one 
quadrant  of  the  computational  domain,  the  actuated  region  within  the  6x6u,  12xl2u  and 
16xl6u  meshes  is  discretized  by  uniform  groups  of  3x3  elements,  6x6  elements  and  8x8 
elements  respectively.  The  two  graded  meshes  include  a  12x12  mesh  (denoted  12xl2g) 
where  the  smallest  element  at  the  actuator  free  edge  has  a  width  of  P/40,  and  a  16x16 
mesh  (denoted  16xl6g)  where  the  smallest  element  (at  the  actuator  edge)  has  a  width  of 
P/80.  These  two  graded  meshes  are  created  in  an  attempt  to  provide  increased  resolution 
at  the  actuator  free  edge  without  increasing  the  problem  size  beyond  that  represented  in 
the  uniform  meshes  12xl2u  and  16x1 6a 

It  should  be  emphasized  that  certain  types  of  physical  phenomena  such  as  damage, 
plasticity  and  wave  propagation  are  often  observed  to  exhibit  non- smooth  distributions 
with  a  high  degree  of  localization.  For  problems  involving  these  types  of  phenomena,  the 
use  of  low  order  finite  elements  (particularly  linear  elements)  is  generally  favored  in 
practice  since  they  minimize  the  domain  over  which  a  localized  behavior  is  dispersed. 
However,  in  the  modeling  of  structural  members  that  exhibit  significant  bending  behavior 
(e.g.  beams,  plates,  shells),  quadratic  elements  are  known  to  be  superior  to  linear 
elements  since  the  latter  exhibits  a  much  higher  degree  of  spurious  transverse  shear 
stiffness27.  Therefore,  despite  the  fact  that  quadratic  elements  are  not  as  appropriate  for 
modeling  localized,  non-smooth  behavior  as  linear  elements,  the  present  study  is 
conducted  using  quadratic  elements  due  to  their  ability  to  deliver  higher  quality  strain 
distributions  for  bending  dominated  problems. 

Regardless  of  the  layerwise  model  type  or  the  level  of  2-D  mesh  density,  all  solutions 
are  computed  using  the  same  modest  level  of  transverse  discretization  that  consists  of  two 
equal  linear  layers  for  each  of  the  piezoceramic  actuators,  one  linear  layer  for  each  of  the 
adhesive  bond  layers  and  four  equal  linear  layers  for  the  aluminum  plate.  Since  the 
adhesive  material  is  modeled  with  a  single  linear  layer,  the  computed  transverse  shear 
strain  and  transverse  normal  strain  in  the  adhesive  layer  necessarily  represent  the 
thickness  average  values  of  these  quantities  in  the  adhesive  layer.  This  particular  level 
of  transverse  discretization  is  chosen  because  it  is  coarse  enough  to  be  practical  for  larger 
problems,  while  still  permitting  a  reasonably  accurate,  piecewise- constant  representation 
of  the  transverse  strains  through  the  thickness  of  the  three  materials. 

Figure  5  shows  a  representative  selection  of  transverse  shear  strain  and  transverse 
normal  strain  distributions  in  the  adhesive  layer  caused  by  imposing  an  actuation  strain 
magnitude  of  £i=e 2=0.001  in  the  piezoceramic  patches.  The  transverse  strains  shown  in 
Figure  5  are  computed  along  a  diagonal  line  that  runs  from  the  center  of  the  actuator  to 
the  free  corner  of  the  actuator.  This  particular  distribution  is  shown  since  the  free  comer 
of  the  adhesive  layer  exhibits  the  most  severe  local  strain  state.  Representative  linear 
elastic  results  are  shown  for  both  the  LW1  and  LW2  models  using  all  five  levels  of  2-D 
mesh  density,  specifically,  the  transverse  shear  strain  distribution  is  shown  for  the  LW 1 
model,  and  the  transverse  normal  strain  distribution  is  shown  for  the  LW2  model.  In 
each  case,  the  strain  values  are  computed  only  at  the  reduced  Gaussian  integration  points 
within  each  element  and  these  computed  values  are  marked  by  various  symbols  in  parts  A 
through  D  of  Figure  5 .  Parts  A  and  B  of  Figure  5  show  results  for  the  thicker  actuated 
plate  configuration  (2P/H=2),  while  parts  C  and  D  show  the  results  for  the  thinner 
actuated  plate  configuration  (2P/H=8). 


21 


Figure  5.  Comparison  of  adhesive  strain  distribution  near  the  free  corner  of  the  upper 
actuator  patch  as  predicted  by  linear  elastic  LW1  and  LW2  models  for  the  case  of 
bending  actuation  of  a  simply  supported  plate. 


As  seen  in  parts  A  through  Dof  Figure  5,  all  five  levels  of  2-D  mesh  density  are 
capable  of  producing  a  local  adhesive  response  with  qualitatively  correct  characteristics. 
In  each  case,  the  transverse  strains  are  insignificantly  small  in  the  interior  region  of  the 
actuator  and  increase  dramatically  as  the  free  corner  of  the  actuator  is  approached.  To  aid 
visualization  of  the  transverse  strain  distributions,  the  individual  computed  values  of  the 
most  refined  mesh  (16xl6g)  are  connected  by  straight  lines.  The  main  difference 
between  the  computed  results  from  the  two  thickness  configurations  (2P/H=2  vs. 

2P/H=8)  is  the  width  of  the  boundary  layer  region  over  which  the  transverse  strains  attain 
significant  magnitude.  For  the  thicker  configuration  (2P/H=2),  the  width  of  the  adhesive 
boundary  layer  region  is  approximately  0.5P,  while  for  the  thinner  configuration 
(2P/H=8),  the  width  of  the  adhesive  boundary  layer  region  is  approximately  0.15P. 

For  the  thicker  configuration  (2P/H=2),  the  results  from  the  12xl2u,  12xl2g,  16xl6u 
and  16xl6g  meshes  show  very  good  agreement  with  all  computed  strain  values  lying  on  a 
smooth  curve,  while  the  results  from  the  6x6u  mesh  can  are  seen  to  exhibit  a  mild 
oscillation  about  the  smooth  curve.  For  the  thinner  configuration  (2P/H=8),  the  width  of 
the  boundary  layer  region  is  too  narrow  to  be  adequately  resolved  by  either  of  the  three 
uniform  meshes  (6x6u,  12xl2u,  16xl6u);  however,  the  results  from  the  two  graded 
meshes  (12xl2g  and  16xl6g)  show  very  good  agreement  with  all  computed  values  lying 
on  a  single  smooth  curve. 


22 


Based  on  the  results  shown  in  Figure  5,  the  12xl2g  mesh  is  observed  to  offer  the  best 
compromise  between  solution  quality  and  computational  efficiency.  It  should  be  noted 
that  the  16xl6g  mesh  consistently  produces  peak  transverse  strains  that  are  slightly 
higher  than  those  produced  by  the  12xl2g  mesh;  however,  this  difference  is  not 
anticipated  to  significantly  affect  the  results  of  a  progressive  damage  analysis,  since  it 
simply  results  in  a  slightly  earlier  initiation  of  local  damage  in  the  16xl6g  mesh.  Once 
local  damage  is  initiated,  the  strain  distributions  no  longer  exhibit  the  sharp  peaks  as  seen 
in  Figure  5,  thus  the  extra  local  resolution  of  the  16xl6g  mesh  is  not  warranted  given  its 
increased  computational  cost. 

Assuming  then  that  the  12xl2g  mesh  is  adequate  for  providing  the  necessary 
resolution  of  the  transverse  strains  in  the  adhesive  boundary  layer  region,  Figure  6  shows 
a  comparison  of  the  linear  elastic  stress  distributions  predicted  by  the  LW1  and  LW2 
models  on  the  12xl2g  mesh.  Examination  of  Figure  6  reveals  that  for  both  thickness 
configurations  (2P/H=2  and  2P/H=8),  the  LW1  model  and  LW2  model  predict  similar 
transverse  shear  stress  distributions  in  the  adhesive  layer,  with  the  LW1  predicting 
slightly  higher  transverse  shear  stress  than  the  LW2  model.  However,  the  LW1  model  is 
based  on  the  assumption  of  zero  transverse  normal  stress,  while  the  LW2  model  predicts 
that  the  peak  transverse  normal  stress  is  approximately  twice  the  magnitude  of  the  peak 
transverse  shear  stress.  Normally,  neglecting  the  transverse  normal  stress  is  an  accepted 
component  of  most  beam,  plate  and  shell  models;  however,  in  the  remainder  of  this 
study,  the  explicit  inclusion  of  transverse  normal  stress  is  shown  to  have  a  profound 
effect  on  local  damage  progression.  Finally,  the  LW1  and  LW2  models  differ  in  the 
predicted  local  distribution  of  inplane  normal  stress  within  the  boundary  layer  region, 
namely,  the  LW2  model  shows  a  marked  increase  in  inplane  normal  stress  as  the  free 
comer  of  the  actuator  is  approached,  while  the  LW 1  model  shows  very  little  change  in 
inplane  normal  stress  as  the  free  comer  of  the  actuator  is  approached 
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5.2  Progressive  damage  of  the  adhesive  layer  under  imposed  bending  actuation 

Now  let  us  consider  the  local  damage  evolution  that  occurs  in  the  adhesive  layer  of 
the  simply- supported  actuated  plate  when  subjected  to  bending  actuation  alone,  without 
any  other  form  of  mechanical  loading.  In  this  case,  the  top  and  bottom  actuators  are 
subjected  to  opposite  voltages  of  equal  magnitude  causing  the  plate  to  bend.  The 
objective  is  to  determine  the  mode,  severity  and  extent  of  damage  in  the  adhesive  bond 
layer  as  influenced  by  the  choice  of  model  type  (LW1  vs.  LW2)  and  the  thickness 
configuration  of  the  actuated  plate.  To  this  end,  two  different  thickness  configurations 
(2P/H=2  and  2P/H=8)  are  simulated  using  both  LW1  and  LW2  models.  Based  on  the 
results  of  the  previous  section,  the  12xl2g  2-D  mesh  is  utilized  in  each  case. 

Furthermore,  the  present  analyses  make  use  of  the  same  level  of  thickness  discretization 
used  in  the  previous  section,  namely,  two  equal  linear  layers  for  each  of  the  piezoceramic 
actuators,  one  linear  layer  for  each  of  the  adhesive  bond  layers  and  four  equal  linear 
layers  for  the  aluminum  plate. 

The  aluminum  material  and  piezoceramic  material  are  modeled  as  linear  elastic 
materials,  while  the  adhesive  material  is  permitted  to  exhibit  progressive  damage 
(material  nonlinearity).  The  damage  surface  constants  and  the  damage  hardening 
constants  that  are  assumed  for  the  adhesive  material  are  given  below  (same  as  used  in 
Section  3.4). 

Adhesive  damage  surface  constants'.  Jn  =  J22  =  J33  =  1, 

Adhesive  damage  hardening'.  y(|3)  =  2.5(1 06)()  and  y(0)  =  0. 

For  both  thickness  configurations,  the  peak  electric  field  applied  to  each  actuator  is 
sufficient  to  produce  inplane  normal  strains  of  8 1=82=0.001  in  an  unconstrained 
piezoceramic  patch  This  peak  electric  field  is  achieved  by  the  application  of  a  series  of 
20  equal  load  increments  and  is  chosen  because  it  represents  the  upper  limit  of  actuation 
strain  that  can  be  obtained  with  a  typical  piezoceramic  material. 

For  the  maximum  imposed  actuation  strain  (0.001),  Parts  A  through  D  of  Figure  7 
show  the  distribution  of  damage  eigenvalues  in  the  adhesive  layer  as  predicted  by  the 
LW1  and  LW2  models  for  both  thickness  configurations.  Parts  A  and  C  show  the 
adhesive  damage  eigenvalues  D|  and  Eh  plotted  along  the  x  axis  (normalized  coordinate 
0  <  r/P  <  1)  for  the  thickness  configurations  2P/H=2  and  2P/H=8  respectively.  Both  the 
LW1  model  and  the  LW2  model  predict  that  significant  levels  of  adhesive  damage  occur 
only  within  the  boundary  layer  region  with  the  maximum  damage  occurring  at  the  free 
edge  of  the  adhesive  (r/P  =1).  Furthermore,  both  models  predict  slightly  higher  damage 
for  the  thicker  configuration  (2P/H=2)  than  for  the  thinner  configuration  (2P/H=8),  and 
both  models  predict  peak  local  damage  values  of  less  than  0.10  along  the  x  axis. 

However,  the  two  models  show  some  disagreement  on  the  relative  magnitudes  of  the 
damage  eigenvalues.  For  example,  the  LW2  model  predicts  Eb>Di  while  the  LW1  model 
predicts  D3~Di  .  This  can  be  explained  by  the  fact  that  the  LW 1  model  does  not  include 
transverse  normal  stress;  therefore,  local  adhesive  damage  is  primarily  driven  by  the 
transverse  shear  stress  oxz  which  contributes  equally  to  D|  and  D3  in  the  adhesive.  In 
contrast,  local  adhesive  damage  in  the  LW2  model  is  driven  by  both  transverse  shear 
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stress  Gxz  and  transverse  normal  stress  ozz.  Since  the  transverse  normal  stress  ozz 
contributes  to  D3  but  not  to  Di,  the  LW2  model  predicts  that  D3>Di. 

Parts  B  and  D  of  Figure  7  show  the  adhesive  damage  eigenvalues  Di,  Do  and  Eh 
plotted  along  the  free  edge  of  the  adhesive  layer  (normalized  coordinate  0  <  e/P  <  1)  for 
the  thickness  configurations  2P/H=2  and  2P/H=8  respectively.  Both  the  LW1  model  and 
the  LW2  model  predict  that  adhesive  damage  is  relatively  constant  over  much  of  the  free 
edge  of  the  adhesive  but  then  increases  dramatically  as  the  free  comer  of  the  adhesive  is 
approached  (i.e.  as  e/P— >1).  This  increased  damage  near  the  free  comer  is  caused  by 
higher  stress  concentration  that  results  from  the  two  free  edges  that  meet  at  a  right  angle. 
At  the  free  corner  of  the  adhesive,  both  the  LW1  and  LW2  models  predict  that  D3  is 
much  larger  than  Di  or  D2,  thus  the  local  damaged  state  clearly  suggests  the  early  stages 
of  a  local  delamination  (debonding),  i.e.,  distributed  microcracks  that  are  oriented 
perpendicular  to  the  transverse  (z)  direction  .  However,  the  value  of  D3  predicted  by  the 
LW2  model  is  significantly  higher  than  the  value  of  D3  predicted  by  the  LW1  model. 
Again,  the  LW2  model’s  higher  value  of  D?  is  caused  by  the  fact  that  the  LW2  model 
explicitly  includes  transverse  normal  stress  which  achieves  a  very  high  tensile  value  near 
the  free  edge  and  makes  a  significant  contribution  to  the  growth  of  Eh-  Both  models 
predict  that  Di=D2  at  the  free  comer  (e/P=l)  which  can  be  understood  by  noting  that  the 
diagonal  line  x=y  is  an  axis  of  symmetry  which  results  in  oxx  =  oyy  and  oxz  =  oyz  at  e/P=l. 
Thus  the  stress  components  that  contribute  to  Di  are  equal  to  the  stress  components  that 
contribute  to  Db. 
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Figure  7.  Comparison  of  damage  distribution  in  the  adhesive  layer  as  predicted  by  LW1 
and  LW2  models  for  the  case  of  bending  actuation  of  a  simply  supported  plate.  Imposed 
actuation  strain  =  0.001. 


26 


Given  the  fact  that  the  highest  adhesive  damage  occurs  at  the  free  corner  of  the 
adhesive  layer,  Figures  7  through  9  show  the  predicted  damage  distributions  and 
predicted  stress  distributions  in  the  adhesive  layer  along  a  diagonal  line  that  runs  from  the 
center  of  the  actuated  region  to  the  free  comer  of  the  actuator  (normalized  coordinate  0  < 
s/P'  <  1).  Parts  A  and  B  of  Figure  8  show  the  distribution  of  damage  eigenvalues  in  the 
thicker  configuration  (2P/H=2)  at  various  levels  of  induced  actuator  strain,  as  predicted 
by  both  the  LW 1  and  LW2  models.  Examination  of  Figure  8  reveals  several  noteworthy 
trends.  First,  both  models  predict  that  damage  is  highest  at  the  free  comer  of  the  adhesive 
layer  (i.e.  at  s/P  -1)  and  becomes  insignificant  as  the  center  of  the  actuator  is  approached 
(i.e.,  as  s/P'— >0).  Second,  both  models  predict  that  the  damage  eigenvalue  D3  is 
significantly  higher  than  damage  eigenvalues  Di  and  D2,  thus  the  primary  form  of 
adhesive  damage  consists  of  microscopic  cracks  that  are  parallel  to  the  reference  surface 
of  the  actuated  plate  and  represent  the  early  stages  of  a  local  delamination  (or  debonding) 
of  the  actuator. 

The  most  striking  difference  between  the  LW1  and  LW2  models  is  the  predicted 
magnitude  of  the  damage  eigenvalues.  The  LW2  model  predicts  that  the  dominant 
damage  eigenvalue  (D3)  is  almost  twice  as  large  as  predicted  by  the  LW  1  model  (see  Part 
A  of  Figure  8).  This  increased  growth  of  D3  in  the  LW2  model  is  a  direct  consequence  of 
explicitly  including  transverse  normal  strain  and  transverse  normal  stress  which  obtain 
very  high  tensile  values  near  the  free  comer  of  the  adhesive.  Consequently,  in  the  LW2 
model,  local  damage  evolution  in  the  adhesive  layer  is  driven  by  both  transverse  shear 
and  transverse  normal  deformation,  while  the  local  adhesive  damage  in  the  LW  1  model  is 
driven  only  by  transverse  shear  deformation.  In  contrast,  the  LW  1  model  predicts  larger 
magnitudes  for  the  non-dominant  damage  eigenvalues  (Di  and  D2)  than  the  LW2  model. 
This  last  observation  is  due  partly  to  the  fact  that  the  LW  1  model  inherently  predicts 
slightly  higher  transverse  shear  deformation  than  the  LW2  model,  and  partly  due  to  the 
fact  that  the  LW2  model’s  increased  damage  eigenvalue  D3  forces  some  of  the  shear 
forces  to  be  redistributed,  thus  limiting  the  growth  of  Di  and  EL 
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Figure  8.  Comparison  of  damage  distribution  in  the  adhesive  layer  as  predicted  by  LW1 
and  LW2  models  for  the  case  of  bending  actuation  of  a  simply  supported  plate  with 
2P/H=2  and  ha//?=0. 1 0.  Results  are  shown  at  four  different  actuation  strain  levels. 
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Figure  9.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  LW1 
and  LW2  models  for  the  case  of  bending  actuation  of  a  simply  supported  plate  with 
2P/H=2  and  hjh=0. 1 0.  Results  are  shown  at  four  different  actuation  strain  levels. 


Figure  9  shows  the  distribution  of  transverse  shear  stress  and  transverse  normal  stress 
in  the  damaged  adhesive  material  of  the  thicker  actuated  plate  (2P/H=2).  In  each  case, 
the  damaged  stresses  are  computed  only  at  the  reduced  Gaussian  integration  points  within 
the  adhesive  material  and  are  plotted  along  a  diagonal  line  that  mns  from  the  center  of  the 
actuated  region  to  the  free  corner  of  the  actuator.  Examination  of  Part  A  of  Figure  9 
reveals  that  the  LW2  model  predicts  lower  values  of  transverse  shear  stress  than  the  LW1 
model,  especially  at  higher  load  levels.  For  low  load  levels  (e.g.  actuation  strain  < 
0.00025),  the  level  of  material  damage  is  minimal;  consequently,  the  difference  between 
the  FW 1  and  FW2  transverse  shear  stress  is  similar  to  that  observed  earlier  for  the  linear 
elastic  case.  However,  at  high  load  levels  (e.g.  actuation  strain  >  0.00075)  the  difference 
in  predicted  transverse  shear  stress  is  quite  noticeable.  This  last  observation  is  due 
primarily  to  the  advanced  state  of  damage  (D3)  predicted  by  the  FW2  model  which 
severely  limits  the  amount  of  transverse  shear  stress  (see  Part  A)  and  transverse  normal 
stress  (see  Part  B)  that  the  material  can  support. 

For  the  thicker  actuated  plate  (2P/H=2),  Figure  10  shows  the  difference  between  the 
adhesive  stresses  obtained  by  linear  elastic  and  progressive  damage  solutions  using  the 
FW1  and  FW2  models.  The  stresses  in  Figure  10  are  computed  at  the  highest  load  level 
(actuation  strain  =  0.001).  The  linear  elastic  solutions  show  significant  stress 
concentrations  at  the  free  comer  of  the  adhesive,  while  the  stresses  obtained  from  the 
progressive  damage  solutions  are  limited  in  magnitude  by  the  locally  weakened  state  of 
the  adhesive  material. 


Figure  10.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  linear 
elastic  and  progressive  damage  solutions  obtained  with  FW1  and  FW2  models.  Results 
shown  for  a  simply  supported  plate  subjected  to  bending  actuation  (actuation  strain  = 
0.001).  2P/H=2and  V/i=0.10. 

Results  for  the  thinner  actuated  plate  (2P/H=8)  are  shown  in  Figures  10  through  12 
and  are  entirely  analogous  to  the  previous  Figures  7  through  9  respectively.  The  damage 
distributions  shown  in  Figure  1 1  and  the  transverse  stress  distributions  shown  in  Figure 
12  show  the  same  basic  trends  for  the  thinner  actuated  plate  as  observed  earlier  for  the 
thicker  actuated  plate.  Figure  13  shows  a  comparison  between  linear  elastic  and 
progressive  damage  solutions  for  the  highest  load  level  (actuation  strain  =  0.001).  In 
comparing  the  results  obtained  for  the  two  thickness  configurations  (2P/H=2  and 
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2P/H=8),  the  main  difference  is  that  the  boundary  layer  region  is  much  narrower  for  the 
thinner  configuration;  consequently,  the  boundary  layer  region  of  the  thinner 
configuration  is  discretized  with  fewer  elements  than  the  thicker  configuration,  thus  the 
damage  and  stress  distributions  seen  in  Figures  10  through  12  do  not  have  the  same 
degree  of  smoothness  as  seen  earlier  in  Figures  7  through  9.  However,  despite  the  fact 
that  fewer  elements  are  used  to  resolve  the  boundary  layer  region  in  the  thinner 
configuration,  the  predicted  behavior  is  similar  for  both  thickness  configurations. 


Figure  11.  Comparison  of  damage  distribution  in  the  adhesive  layer  as  predicted  by 
LW 1  and  LW2  models  for  the  case  of  bending  actuation  of  a  simply  supported  plate  with 
2P/H=8  and  hjh=0. 1 0.  Results  are  shown  at  four  different  actuation  strain  levels. 


Figure  12.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  LW 1 
and  LW2  models  for  the  case  of  bending  actuation  of  a  simply  supported  plate  with 
2P/H=2  and  hjh=0. 1 0.  Results  are  shown  at  four  different  actuation  strain  levels. 
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Figure  13.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  linear 
elastic  and  progressive  damage  solutions  obtained  with  LW1  and  LW2  models.  Results 
shown  for  a  simply  supported  plate  subjected  to  bending  actuation  (actuation  strain  = 
0.001).  2P/H=8  and  \lh= 0.10. 

Finally,  let  us  consider  the  effect  of  local  adhesive  damage  on  the  ability  of  the 
actuators  to  produce  global  deformation  in  the  actuated  plate.  Table  1  shows  the 
predicted  normalized  center  deflection  of  the  simply  supported  plate  under  bending 
actuation  (imposed  actuation  strain  =  0.001).  The  center  deflections  are  obtained  from 
both  linear  elastic  solutions  and  progressive  damage  solutions  using  both  layerwise 
model  types  and  both  thickness  configurations.  Examination  of  the  values  in  Table  1 
reveals  that,  for  all  cases  considered,  the  localized  adhesive  damage  reduces  the  center 
deflection  by  less  than  1%.  Thus  it  is  concluded  that,  although  the  adhesive  damage  is 
locally  severe,  it  is  not  pervasive  enough  to  seriously  degrade  the  overall  effectiveness  of 
the  actuators. 


model 

type 

2P/H 

imposed 

actuation 

strain 

solution 

type 

center 

deflection 

w(0,0,0)/hp 

%  difference  in 
w(0,0,0) 

(damage  vs.  elastic) 

LW1 

2 

0.001 

progressive  damage 

-0.002485 

LW1 

2 

0.001 

linear  elastic 

-0.002502 

0.67 

LW2 

2 

0.001 

progressive  damage 

-0.002431 

LW2 

2 

0.001 

linear  elastic 

-0.002450 

0.80 

LW1 

8 

0.001 

progressive  damage 

-0.045378 

LW1 

8 

0.001 

linear  elastic 

-0.045436 

0.13 

LW2 

8 

0.001 

progressive  damage 

-0.044967 

LW2 

8 

0.001 

linear  elastic 

-0.045035 

0.15 

Table  1.  Comparison  of  center  deflection  as  influenced  by  model  type,  solution  type  and 
thickness  configuration. 

5.3  Progressive  adhesive  damage  (bending  actuation  plus  transverse  loading) 

The  previous  example  considered  dama  ge  to  the  adhesive  bond  layers  that  resulted 
from  bending  actuation  alone,  without  any  other  form  of  mechanical  loading.  It  was 
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demonstrated  that  an  imposed  actuation  strain  of  0.001  was  sufficient  to  produce 
significant  damage  in  the  adhesive  bond,  but  this  damage  was  limited  to  the  boundary 
layer  region  where  very  high  values  of  transverse  shear  stress  and  transverse  normal 
stress  occur.  It  was  further  demonstrated  that  although  the  predicted  adhesive  damage 
was  locally  quite  severe,  the  adhesive  damage  was  not  pervasive  enough  to  significantly 
reduce  the  effectiveness  of  the  actuators.  In  this  section,  we  consider  a  more  severe  two- 
stage  loading  of  the  simply  supported  plate.  More  precisely,  the  plate  is  first  loaded  by  a 
high  level  of  bending  actuation.  Then,  while  holding  the  bending  actuation  constant,  the 
plate  is  subjected  to  a  uniform  distributed  transverse  load  that  opposes  the  bending 
actuation.  The  magnitude  of  the  uniform  distributed  load  is  incrementally  increased  until 
it  effectively  cancels  the  global  deformation  caused  by  the  bending  actuation.  This 
results  in  a  plate  that  exhibits  minimal  global  deformation,  but  high  stress  levels. 

Two  different  thickness  configurations  are  considered  (2P/H=2  and  2P/H=8).  In  each 
case,  the  adhesive  layer  thickness  is  10%  of  the  actuator  thickness  (h//z=0.1),  while  the 
actuator  thickness  is  25%  of  the  aluminum  plate  thickness  (h! H=0.25).  In  each  case,  the 
constant  bending  actuation  is  brought  about  by  applying  an  electric  field  that  is  sufficient 
to  cause  inplane  normal  strains  of  e  1=82=0.001  in  an  unconstrained  piezoceramic  patch. 
Note  that  this  is  equivalent  to  the  maximum  level  of  bending  actuation  imposed  in  the 
previous  example ;  therefore,  solutions  for  the  present  problem  are  computed  by  're¬ 
starting'  the  previous  solutions  with  the  addition  of  a  uniform  distributed  transverse  load. 
For  the  thicker  configuration  (2P/H=2),  the  magnitude  of  the  uniform  distributed 
transverse  load  is  increased  in  increments  of  1  MPa  until  the  loading  is  sufficient  to 
cancel  the  bending  actuation.  For  the  thinner  configuration  (2P/H=8),  the  magnitude  of 
the  uniform  distributed  transverse  load  is  increased  in  increments  of  0.0625  MPa  until  the 
loading  is  sufficient  to  cancel  the  bending  actuation.  For  both  cases,  solutions  are 
computed  using  both  the  LW1(2L/1L/4L)  and  LW2(2L/1L/4L)  models  on  the  12xl2g  2- 
D  mesh. 

Figure  14  shows  the  center  deflection  of  the  simply  supported  plate  at  each  loadstep, 
where  loadstep  0  corresponds  to  the  application  of  bending  actuation  (imposed  actuation 
strain  =  0.001)  without  any  distributed  transverse  load.  In  each  subsequent  loadstep,  the 
opposing  uniform  distributed  transverse  load  is  increased  by  1  MPa  for  the  thicker 
configuration  (2P/H=2),  or  0.0625  MPa  for  the  thinner  configuration  (2P/H=8).  For 
comparison,  both  linear  elastic  solutions  and  progressive  damage  solutions  are  obtained 
using  the  LW1  model  and  the  LW2  model.  Part  A  of  Figure  14  shows  that,  for  the 
thicker  configuration,  a  uniform  distributed  transverse  load  of  approximately  4  MPa  is 
required  to  cancel  the  global  bending  caused  by  the  initial  bending  actuation. 

Furthermore,  it  is  noticed  that  the  overall  load- deflection  curves  appear  to  be  linear  for 
both  the  LW1  and  LW2  models  despite  the  fact  that  local  adhesive  damage  continues  to 
accumulate  with  each  loadstep.  However,  the  slope  of  the  load- deflection  curve  is 
shallower  for  the  LW2  than  for  the  LW1  model.  For  the  thinner  configuration.  Part  B  of 
Figure  14  shows  that  a  uniform  distributed  transverse  load  of  approximately  0.375  MPa  is 
required  to  cancel  the  global  bending  caused  by  the  initial  bending  actuation.  In  this 
case,  the  load- deflection  curves  of  the  LW1  and  LW2  models  are  indistinguishable  from 
each  other. 
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Figure  14.  Normalized  center  deflection  produced  by  constant  bending  actuation 
(actuation  strain  =  0.001)  followed  by  an  opposing  uniform  distributed  transverse  load. 


For  the  thicker  actuated  plate  configuration  (2P/H=2),  Figure  15  shows  the  evolution 
of  local  adhesive  damage  that  occurs  as  the  distributed  load  is  increased  from  0  MPa  to  6 
MPa.  Specifically,  Parts  A  and  B  of  Figure  15  show  the  evolution  of  damage  eigenvalues 
D3  and  D1=D2  respectively  as  predicted  by  the  LW1  model.  Parts  C  and  D  of  Figure  15 
show  the  evolution  of  the  damage  eigenvalues  D3  and  Di=D2  respectively  as  predicted  by 
the  LW2  model.  Examination  of  Figure  15  reveals  several  noteworthy  trends.  First,  both 
models  predict  that  the  dominant  form  of  damage  consists  of  distributed  microcracks  that 
are  oriented  perpendicular  to  the  transverse  direction  (characterized  by  Eh);  however,  the 
LW2  model  predicts  significantly  higher  values  of  D3  than  the  LW1  model.  In  fact,  if  the 
distributed  load  is  increased  beyond  6  MPa  (i.e.  beyond  loadstep  6),  the  LW2  model 
predicts  a  local  delamination  failure  (D3=l)  at  the  reduced  Gauss  point  closest  to  the  free 
corner  of  the  adhesive.  Again,  the  LW2  model’s  larger  values  of  D3  can  be  explained  by 
the  fact  that  the  LW2  model  explicitly  includes  transverse  normal  stress  which  is  a 
significant  contributor  to  D3  in  this  particular  problem.  Second,  despite  the  fact  that  the 
LW2  model  predicts  higher  values  of  D3,  both  models  predict  approximately  the  same 
growth  rate  of  D3  as  the  distributed  load  is  incremented.  For  example,  over  the  course  of 
six  distributed  load  increments,  both  models  show  that  the  peak  value  of  E>?  is 
approximately  doubled.  Third,  the  LW 1  model  predicts  larger  peak  values  of  the  non¬ 
dominant  damage  eigenvalues  (Di  and  Eh)  than  the  LW2  model;  moreover,  the  two 
models  predict  qualitatively  different  evolution  characteristics  for  these  non-dominant 
damage  eigenvalues.  The  LW1  model  predicts  that  Di  and  D2  continue  to  exhibit  an 
ever-increasing  maxima  at  the  free  corner  of  the  adhesive.  In  contrast,  the  LW2  model 
predicts  that  the  Dj  and  Eh  maxima  begin  to  move  inward,  away  from  the  free  comer  of 
the  adhesive  as  D3  continues  to  increase. 

Figure  16  shows  the  evolution  of  the  local  transverse  stress  components  in  the 
adhesive  layer  of  the  thicker  configuration  (2P/H=2)  as  the  distributed  load  is  increased 
from  0  MPa  to  6  MPa.  Specifically,  Part  A  of  Figure  16  shows  the  evolution  of  the 
transverse  shear  stress  components  Gxz=ayz  as  predicted  by  the  LW1  model.  Parts  B  and 
C  of  Figure  16  show  the  evolution  of  the  transverse  shear  stress  components  and  the 
transverse  normal  stress  component  respectively  as  predicted  by  the  LW2  model. 
Examination  of  Parts  A-C  of  Figure  16  reveals  that  as  local  damage  accumulates  near  the 
free  corner  of  the  adhesive,  local  load  redistribution  occurs;  i.e.,  the  portion  of  the  stress 
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that  can  not  be  supported  by  the  damaged  material  is  transferred  to  neighboring  material 
that  is  less  damaged.  In  Parts  A-C  of  Figure  16,  as  the  distributed  load  is  increased,  the 
location  of  the  transverse  stress  maxima  shifts  to  the  left  into  lesser  damaged  material.  At 
the  same  time,  the  stress  near  the  free  comer  of  the  adhesive  decreases  commensurate 
with  the  increasing  level  of  local  damage.  Also,  it  is  noted  that  within  any  particular  plot 
(A,  B,  or  C),  the  peak  stress  value  remains  constant  as  the  stress  distribution  shifts  to  the 
left  (at  least  to  the  extent  that  the  chosen  mesh  can  resolve  the  stress  field).  This  constant 
peak  stress  is  consistent  with  a  material  that  can  only  support  a  limited  level  of  stress  for 
any  given  deformation  state.  Interestingly,  the  LW2  model  predicts  that  the  adhesive 
material  can  support  a  maximum  transverse  shear  stress  of  cxz=cyz=22.5  MPa,  while  the 
LW 1  model  predicts  that  the  adhesive  material  can  support  a  maximum  transverse  shear 
stress  of  oxz=ayz=21 .5  MPa.  It  should  be  noted  that  the  exact  same  damage  evolution 
equations  are  used  in  both  models;  thus,  the  difference  in  peak  shear  stress  is  caused  by 
the  fact  that  the  LW2  model  predicts  more  of  a  local  triaxial  state  of  stress  than  the  LW1. 
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Figure  15.  Comparison  of  damage  distribution  in  the  adhesive  layer  as  predicted  by 
LW1  model  (Parts  A  &  B)  and  LW2  model  (Parts  C  &  D).  Uniform  distributed 
transverse  load  increments  (1  MPa  each)  are  applied  to  oppose  pre-existing  bending 
actuation  (actuation  strain  =  0.001). 
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Figure  16.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  LW 1 
and  LW2  models.  Uniform  distributed  transverse  load  increments  (1  MPa)  are  applied  to 
oppose  pre-existing  bending  actuation  (actuation  strain  =  0.001).  2P/H=2,  \\Jh=(). 1 0. 


For  the  thicker  actuated  plate  configuration  (2P/H=2),  Figure  17  shows  the 
distribution  of  transverse  stresses  in  the  adhesive  layer  at  a  distributed  load  of  4  MPa. 
This  particular  level  of  transverse  loading  is  just  sufficient  to  counteract  the  imposed 
bending  actuation.  The  stresses  in  Figure  17  are  obtained  from  progressive  damage 
solutions  and  linear  elastic  solutions  using  both  the  LW1  and  LW2  models.  As  seen  in 
Figure  17,  the  linear  elastic  solutions  predict  significant  transverse  stress  concentrations 
near  the  free  comer  of  the  adhesive.  In  the  progressive  damage  solutions,  these 
transverse  stress  concentrations  can  not  be  tolerated  and  result  in  adhesive  damage  which 
severely  limits  the  magnitude  of  the  local  transverse  stresses. 
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Figure  17.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  linear 
elastic  and  progressive  damage  solutions  obtained  with  LW1  and  LW2  models.  Results 
shown  for  a  simply  supported  plate  subjected  to  combined  bending  actuation  (actuation 
strain  =  0.001)  and  opposing  uniform  distributed  transverse  load  (4  MPa).  2P/H=2  and 
ha//z=0.10. 

Results  for  the  thinner  actuated  plate  configuration  (2P/H=8)  are  shown  in  Figures  17 
through  19  which  are  entirely  analogous  to  earlier  Figures  14  through  16  respectively. 
Figures  17  and  18  show  the  distribution  of  damage  eigenvalues  and  damaged  transverse 
stresses  predicted  by  the  LW1  and  LW2  models  for  the  thinner  configuration  (2P/H=8) 
and  show  the  same  general  trends  discussed  earlier  for  the  thicker  configuration.  Figure 
20  shows  a  comparison  between  linear  elastic  and  progressive  damage  solutions  at  a 
uniform  distributed  transverse  load  of  0.375  MPa.  This  particular  level  of  transverse 
loading  is  just  sufficient  to  counteract  the  imposed  bending  actuation  (actuation  strain  = 
0.001).  In  comparing  the  results  obtained  for  the  two  thickness  configurations(2P/H=2 
and  2P/H=8),  the  main  difference  is  that  the  boundary  layer  region  is  much  narrower  for 
the  thinner  configuration;  consequently,  the  boundary  layer  region  of  the  thinner 
configuration  is  discretized  with  fewer  elements  than  the  thicker  configuration,  thus  the 
damage  and  stress  distributions  seen  in  Figures  17  through  19  do  not  have  the  same 
degree  of  smoothness  as  seen  earlier  in  Figures  14  through  16.  However,  despite  the  fact 
that  fewer  elements  are  used  to  resolve  the  boundary  layer  region  in  the  thinner 
configuration,  the  predicted  behavior  is  similar  for  both  thickness  configurations. 
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Figure  18.  Comparison  of  damage  distribution  in  the  adhesive  layer  as  predicted  by 
LW1  and  LW2  models.  Uniform  distributed  transverse  load  increments  (0.0625  MPa) 
are  applied  to  oppose  pre-existing  bending  actuation  (actuation  strain  =  0.001).  2P/H=8, 
ha//r=0.10. 
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Figure  19.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  LW 1 
and  LW2  models.  Uniform  distributed  transverse  load  increments  (0.0625  MPa)  are 
applied  to  oppose  pre-existing  bending  actuation  (actuation  strain  =  0.001).  2P/H=8, 
ha//?=0. 1 0. 


Figure  20.  Comparison  of  stress  distribution  in  the  adhesive  layer  as  predicted  by  linear 
elastic  and  progressive  damage  solutions  obtained  with  LW1  and  LW2  models.  Results 
shown  for  a  simply  supported  plate  subjected  to  combined  bending  actuation  (actuation 
strain  =  0.001)  and  opposing  uniform  distributed  transverse  load  (0.375  MPa).  2P/H=8 
and  \\/h=0. 1 0. 
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5.4  Progressive  adhesive  damage  ( distributed  transverse  loading  only ) 

As  a  final  example,  let  us  consider  the  evolution  of  local  adhesive  damage  that  occurs 
as  the  actuated  plate  is  subjected  to  a  uniform  distributed  transverse  load  without  any 
electrical  field  applied  to  the  actuators.  Again,  two  different  thickness  configurations  are 
considered  (2P/H=2  and  2P/H=8).  In  each  case,  the  adhesive  layer  thickness  is  10%  of 
the  actuator  thickness  (iy/z=0.1),  while  the  actuator  thickness  is  25%  of  the  aluminum 
plate  thickness  (hi H=0.25).  In  both  cases,  the  uniform  distributed  transverse  load  is 
applied  to  the  bottom  surface  of  the  actuated  plate  and  acts  in  an  upward  direction.  For 
the  thicker  configuration  (2P/H=2),  the  magnitude  of  the  uniform  distributed  transverse 
load  is  increased  in  increments  of  1  MPa  until  the  loading  is  sufficient  to  initiate  a 
localized  material  failure  in  one  of  the  adhesive  layers.  For  the  thinner  configuration 
(2P/H=8),  the  magnitude  of  the  uniform  distributed  transverse  load  is  increased  in 
increments  of  0.0625  MPa  until  the  loading  is  sufficient  to  initiate  a  localized  material 
failure  in  one  of  the  adhesive  layers.  For  both  cases,  solutions  are  computed  using  the 
LW2(2L/1L/4L)  model  on  the  12xl2g  2-D  mesh. 

For  the  thicker  configuration  (2P/H=2),  the  initial  material  failure  occurred  at  a  load 
of  15  MPa,  while  the  thinner  configuration  exhibited  the  initial  material  failure  at  a  load 
of  1.3125  MPa.  In  both  cases,  the  initial  material  failure  occurred  at  the  reduced  Gauss 
point  located  closest  to  the  free  comer  of  the  lower  adhesive  layer  and  was  indicative  of  a 
delamination  initiation,  or  debonding  (i.e.  E>3=1).  Figure  21  shows  the  distribution  of  the 
damage  eigenvalues  and  the  distribution  of  the  transverse  stresses  that  occur  in  the  lower 
adhesive  layers  of  both  configurations  during  the  loadstep  immediately  preceding  the 
initial  localized  adhesive  failure.  Examination  of  Parts  A  and  C  of  Figure  21  reveals  that 
the  adhesive  damage  is  primarily  characterized  by  a  high  concentration  of  microcracks 
that  are  oriented  perpendicular  to  the  thickness  direction  (i.e.,  E>3»Di=D2).  Both 
thickness  configurations  exhibit  damage  distributions  with  similar  shapes  and  similar 
ratios  between  the  various  damage  eigenvalues.  The  main  difference  is  simply  the  width 
of  the  boundary  layer  where  significant  adhesive  damage  occurs.  Parts  B  and  D  of  Figure 
21  show  the  distribution  of  the  transverse  stresses  that  occur  in  the  lower  adhesive  layers 
of  both  configurations.  In  order  to  emphasize  the  local  load  redistribution  that  occurs  as 
damage  accumulates  near  the  free  corner  of  the  adhesive,  Parts  B  and  D  of  Figure  21  also 
show  the  distribution  of  the  transverse  stresses  predicted  by  linear  elastic  (non- damaging) 
solutions.  A  comparison  of  these  stress  distributions  shows  that  the  linear  elastic  stress 
concentrations  are  completely  absent  in  the  progressive  damage  solution. 
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Figure  21.  Damage  and  transverse  stress  distribution  in  the  lower  adhesive  bond  layer  for 
simply  supported  actuated  plate  subjected  to  uniform  distributed  transverse  load  with 
passive  actuators. 


Finally,  let  us  consider  the  effect  of  local  adhesive  damage  on  the  overall  structural 
stiffness  of  the  actuated  plate.  In  the  present  example  problem,  the  actuators  are  passive; 
therefore,  the  actuators  simply  provide  an  additional  bending  stiffness  to  the  aluminum 
plate.  However,  if  the  adhesive  bond  layers  are  damaged,  then  the  aluminum  plate  is 
effectively  prevented  from  utilizing  all  of  the  extra  stiffness  provided  by  the  actuators; 
consequently,  the  actuated  plate  exhibits  increased  overall  compliance.  To  quantify  this 
effect,  Table  2  shows  the  transverse  center  deflection  w(0,0,0)  predicted  for  both 
configurations  by  linear  elastic  solutions  and  progressive  damage  solutions.  As  seen  in 
Table  2,  the  center  deflections  are  larger  for  the  progressive  damage  solutions  than  for  the 
linear  elastic  solutions;  however,  the  increase  in  overall  bending  compliance  is  quite 
small.  For  the  thicker  configuration,  the  damaged  adhesive  layer  allows  a  1.4%  increase 
in  center  deflection.  For  the  thinner  configuration,  the  damaged  adhesive  layer  only 
allows  a  0.7%  increase  in  center  deflection.  Thus,  while  the  adhesive  damage  is  locally 
quite  severe  (see  Parts  A  and  C  of  Figure  21),  the  damage  is  not  pervasive  enough  to 
significantly  degrade  the  stiffness  support  function  of  the  actuators. 
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model 

type 

2P/H 

trans. 

load 

(MPa) 

solution 

type 

center 

deflection 

w(0,0,0)/hp 

%  difference  in 
w(0,0,0) 

(damage  vs.  elastic) 

LW2 

2 

14 

linear  elastic 

0.008252 

LW2 

2 

14 

prog,  damage 

0.008369 

1.4% 

LW2 

8 

1.25 

linear  elastic 

0.15061 

LW2 

8 

1.25 

prog,  damage 

0.15165 

0.7% 

Table  2.  Comparison  of  center  deflections  obtained  by  linear  elastic  solutions  and 
progressive  damage  solutions. 

6.  Summary  and  Conclusions 

This  paper  presented  a  3-D  continuum  damage  mechanics  formulation  and  described 
its  numerical  implementation  into  a  hierarchical,  variable- kinematic  finite  element  model 
that  is  specifically  developed  for  composite  laminates.  The  damage  eigenvalues 
represent  the  effective  fractional  reduction  in  load  carrying  area  on  planes  that  are 
perpendicular  to  the  principal  material  directions  (i.e.,  orthotropic  damage)  and  are 
permitted  to  exhibit  a  spatial  variation  within  each  material  ply  that  is  commensurate  with 
the  predicted  ply- level  strain  fields.  The  resulting  model  was  used  to  simulate  the 
accumulation  of  microscopic  damage  in  the  adhesive  bond  layers  of  a  simple  actuated 
plate  (up  to  the  point  of  initial  material  failure  in  the  adhesive).  Solutions  were  obtained 
using  both  type- 1  and  type- II  layerwise  finite  elements  for  simply- supported  actuated 
plates  subjected  to  bending  actuation  and/or  distributed  transverse  loading. 

The  results  obtained  for  the  actuated  plate  problem  clearly  show  that,  even  under  very 
mild  levels  of  external  loading,  the  high  concentration  of  transverse  shear  strain  and 
transverse  normal  strain  that  exists  near  the  free  edges  of  the  adhesive  layer  leads  to  the 
initiation  of  localized  adhesive  damage.  This  damage  manifests  itself  primarily  in  the 
form  of  distributed  microcracks  that  are  oriented  normal  to  the  transverse  direction  and 
was  shown  to  be  most  severe  at  the  free  comers  of  the  adhesive  layer.  As  the  external 
loading  is  increased,  the  localized  damage  within  the  adhesive  boundary  layer  becomes 
quite  severe  and  can  easily  lead  to  local  material  failure;  however,  due  to  the  limited  size 
of  the  damaged  region,  the  degradation  of  overall  actuator  function  was  shown  to  be 
minimal.  In  addition,  the  damage  to  the  adhesive  layer  did  not  significantly  impair  the 
stiffness  support  function  of  the  actuators  when  operating  in  passive  mode. 

A  comparison  of  solutions  obtained  from  type-I  and  type- II  layerwise  models  showed 
that  the  LW2  model  predicted  significantly  higher  values  of  local  adhesive  damage  than 
the  LW 1  model.  Therefore,  it  can  be  concluded  that  transverse  normal  strain  is  a 
significant  contributor  to  local  damage  accumulation  in  adhesive  layer.  Thus  the 
prediction  of  damage  accumulation  in  the  adhesive  layer  is  shown  to  require  the  use  of  a 
full  3-D  macroscopic  model  that  explicitly  accounts  for  all  six  strain  components  on  a 
discrete  layer  basis. 
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